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We review the electronic properties of bilayer graphene, beginning with a description of the tight- 
binding model of bilayer graphene and the derivation of the effective Hamiltonian describing massive 
chiral quasiparticles in two parabolic bands at low energy. We take into account five tight-binding 
parameters of the Slonczewski-Weiss-McClure model of bulk graphite plus intra- and interlayer 
asymmetry between atomic sites which induce band gaps in the low-energy spectrum. The Hartree 
model of screening and band-gap opening due to interlayer asymmetry in the presence of external 
gates is presented. The tight-binding model is used to describe optical and transport properties 
including the integer quantum Hall effect, and we also discuss orbital magnetism, phonons and the 
influence of strain on electronic properties. We conclude with an overview of electronic interaction 
effects. 
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I. INTRODUCTION 

The production by mechanical exfoliation of isolated 
flakes of graphene with excellent conducting properties 
[l[ was soon followed by the observation of an unusual 
sequence of plateaus in the integer quantum Hall effect in 
monolayer graphene 0, Q . This confirmed the fact that 
charge carriers in monolayer graphene are massless chiral 
quasiparticles with a linear dispersion, as described by a 
Dirac-like effective Hamiltonian and it prompted 

an explosion of interest in the field p} . 
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The integer quantum Hall effect in bilayer graphene [8[ 
is arguably even more unusual than in monolayer because 
it indicates the presence of massive chiral quasiparticles 
@ with a parabolic dispersion at low energy. The effec- 
tive Hamiltonian of bilayer graphene may be viewed as 
a generalisation of the Dirac-like Hamiltonian of mono- 
layer graphene and the second (after the monolayer) in a 
family of chiral Hamiltonians that appear at low energy 
in ABC-stacked (rhombohedral) multilayer graphene [2- 
Il5l ]. In addition to interesting underlying physics, bilayer 
graphene holds potential for electronics applications, not 
least because of the possibility to control both carrier 
density and energy band gap through doping or gating 

[IE IBS]. 

Not surprisingly, many of the properties of bilayer 
graphene are similar to those in monolayer 0, [HJ . These 
include excellent electrical conductivity with room tem- 
perature mobility of up to 40, 000 cm 2 V -1 s _1 in air [22j : 
the possibility to tune electrical properties by chan gin g 
the carrier density through gating or doping 0, fl6| ; 
high thermal conductivity with room temperature ther- 
mal conductivity of about 2, 800 Wm" 1 K" 1 [HHj]; me- 
chanical stiffness, strength and flexibility (Young's mod- 
ulus is estimated to be about 0.8 TPa [2|| [26j]); trans- 
parency with transmittance of white light of about 95 % 
p7j : impermeability to gases [Hj]; and the ability to be 
chemically functionalised [29| . Thus, as with monolayer 
graphene, bilayer graphene has potential for future appli- 
cations in many areas [2l| including transparent, flexible 
electrodes for touch screen displays pjOT^high- frequency 
transistors [Hj]; thermoelectric devices |32l|; photonic de- 
vices including plasmonic devices [H[ and photodetectors 
[34[ ; energy applications including batteries [35, 36]; and 
composite materials [13, HH ■ 

It should be stressed, however, that bilayer graphene 
has features that make it distinct from monolayer. The 
low-energy band structure, described in detail in Sec- 
tion [TT1 is different. Like monolayer, intrinsic bilayer has 
no band gap between its conduction and valence bands, 
but the low-energy dispersion is quadratic (rather than 
linear as in monolayer) with massive chiral quasiparti- 
cles [1, Q rather than massless ones. As there arc two 
layers, bilayer graphene represents the thinnest possible 
limit of an intercalated material [ID, H(| . It is possible 
to address each layer separately leading to entirely new 
functionalities in bilayer graphene including the possibil- 
ity to control an energy band gap of up to about 300 meV 
through doping or gating [9l. TloL [l6l - l2"oj . Recently, this 
band gap has been used to create devices - constrictions 
and dots - by electrostatic confinement with gates [39j . 
Bilayer or multilayer graphene devices may also be prefer- 
able to monolayer ones when there is a need to use more 
material for increased electrical or thermal conduction, 
strength [13, HI] , or optical signature [HJ . 

In the following we review the electronic properties 
of bilayer graphene. Section |TT] is an overview of the 
electronic tight-binding Hamiltonian and resulting band 
structure describing the low-energy chiral Hamiltonian 



and taking into account different parameters that cou- 
ple atomic orbitals as well as external factors that may 
change the electron bands by, for example, opening a 
band gap. We include the Landau level spectrum in the 
presence of a perpendicular magnetic field and the corre- 
sponding integer quantum Hall effect. In section IIIII we 
consider the opening of a band gap due to doping or gat- 
ing and present a simple analytical model that describes 
the density-dependence of the band gap by taking into 
account screening by electrons on the bilayer device. The 
tight-binding model is used to describe transport proper- 
ties, section llVl and optical properties, sectionlVl We also 
discuss orbital magnetism in section fVTI phonons and the 
influence of strain in section I VIII Section IVIIII concludes 
with an overview of electronic-interaction effects. Note 
that this review considers Bernal-stacked (also known as 
AB-stacked) bilayer graphene; we do not consider other 
stacking types such as AA-stacked graphene [4(| , twisted 
graphene [4fl446j or two graphene sheets separated by a 
dielectric with, possibly, electronic interactions between 
them E7H5I. 



II. ELECTRONIC BAND STRUCTURE 
A. The crystal structure and the Brillouin zone 

Bilayer graphene consists of two coupled monolayers 
of carbon atoms, each with a honeycomb crystal struc- 
ture. Figure Q] shows the crystal structure of monolayer 
graphene, figure [3] shows bilayer graphene. In both cases, 
primitive lattice vectors ai and &2 may be defined as 



ai= 2' "2" 



a 2 



%/3q 



(f) 



where a = |ai| = |a.2 1 is the lattice constant, the distance 
between adjacent unit cells, a = 2.46 A [53|, Note that 
the lattice constant is distinct from the carbon-carbon 
bond length acc — a/VS = 1.42 A, which is the distance 
between adjacent carbon atoms. 

In monolayer graphene, each unit cell contains two 
carbon atoms, labelled A and B, figure Q^a) • The po- 
sitions of A and B atoms are not equivalent because it 
is not possible to connect them with a lattice vector of 
the form R = niai + rt2&2, where n\ and rii are inte- 
gers. Bilayer graphene consists of two coupled monolay- 
ers, with four atoms in the unit cell, labelled Al, Bl 
on the lower layer and A2, B2 on the upper layer. The 
layers are arranged so that one of the atoms from the 
lower layer Bl is directly below an atom, A2, from the 
upper layer. We refer to these two atomic sites as 'dimer' 
sites because the electronic orbitals on them are cou- 
pled together by a relatively strong interlayer coupling. 
The other two atoms, Al and B2, don't have a counter- 
part on the other layer that is directly above or below 
them, and are referred to as 'non-dimer' sites. Note that 
some authors [lO, 54 56] employ different definitions of 
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FIG. 1. (a) Crystal structure of monolayer graphene with 
A (B) atoms shown as white (black) circles. The shaded 
rhombus is the conventional unit cell, ai and a2 are primitive 
lattice vectors, (b) Reciprocal lattice of monolayer and bilayer 
graphene with lattice points indicated as crosses, bi and Y>2 
are primitive reciprocal lattice vectors. The shaded hexagon 
is the first Brillouin zone with T indicating the centre, and 
K+ and K- showing two non-equivalent corners. 



A and B sites as used here. The point group of the bi- 
layer crystal structure is D 3 d [l^, H3, HH consisting of 
elements ({E, 2C3, 36*2, i, 2Sq, Sad}), and it may be re- 
garded as a direct product of group D 3 ({E, 2C3, 3C 2 }) 
with the inversion group Ci ({E,i}). Thus, the lattice 
is symmetric with respect to spatial inversion symmetry 
(x,y,z) {-x,-y,-z). 

Primitive reciprocal lattice vectors t>i and b 2 of mono- 
layer and bilayer graphene, where ai ■ t>i = a 2 ■ b 2 = 2n 
and ai ■ b 2 = a 2 ■ bi = 0, are given by 



/2tt 2tt \ , /2tt 2tt \ 

As shown in figurc[TIb), the reciprocal lattice is an hexag- 
onal Bravais lattice, and the first Brillouin zone is an 
hexagon. 




FIG. 2. (a) Plan and (b) side view of the crystal structure of 
bilayer graphene. Atoms AI and BI on the lower layer are 
shown as white and black circles, A2, B2 on the upper layer 
are black and grey, respectively. The shaded rhombus in (a) 
indicates the conventional unit cell. 

B. The tight-binding model 

1. An arbitrary crystal structure 

In the following, we will describe the tight-binding 
model [H, [H, Ho| and its application to bilayer graphene. 
We begin by considering an arbitrary crystal with trans- 
lational invariance and M atomic orbitals <f> m per unit 
cell, labelled by index m = 1 . . . M. Bloch states 
$ m (k, r) for a given position vector r and wave vector 
k may be written as 

1 N 

$ m (k, t) = -=Y, e ik - Rm '^m (r - Rn.,0 , (3) 

where N is the number of unit cells, i = 1 . . . N labels 
the unit cell, and 'R m ,i is the position vector of the mth 
orbital in the zth unit cell. 

The electronic wave function tyj (k, r) may be ex- 
pressed as a linear superposition of Bloch states 

M 

* i (k,r)=53^, ni (k)* ro (k,r) ) (4) 

m— 1 

where ipj m are expansion coefficients. There are M dif- 
ferent energy bands, and the energy Ej(k) of the jth 
band is given by Ej(k) = jlH]^ j) / j\^f j) where H 
is the Hamiltonian. Minimising the energy Ej with re- 
spect to the expansion coefficients ipj, m [H, [60] leads to 

Hipj = Ej Sipj , (5) 
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where ijjj is a column vector, ipj = (ipji, ipj2, ■ ■ ■ , ipjAi)- 
The transfer integral matrix H and overlap integral ma- 
trix S are M x M matrices with matrix elements defined 
by 



= ($ 1$ / 

\ "777. ^ m.' 



(6) 



The band energies Ej may be determined from the gen- 
eralised eigenvalue equation ([5]) by solving the secular 
equation 



det (H-EjS) = 0, 



(7) 



where 'det' stands for the determinant of the matrix. 

In order to model a given system in terms of the gener- 
alised eigenvalue problem ([5]) , it is necessary to determine 
the matrices H and S. We will proceed by considering 
the relatively simple case of monolayer graphene, before 
generalising the approach to bilayers. In the following 
sections, we will omit the subscript j on ipj and Ej in 
equation ([5]), remembering that the number of solutions 
is M, the number of orbitals per unit cell. 



2. Monolayer graphene 

Here, we will outline how to apply the tight-binding 
model to graphene, and refer the reader to tutorial-style 
reviews [53, [6(| for further details. We take into account 
one 2p z orbital per atomic site and, as there are two 
atoms in the unit cell of monolayer graphene, figure (Ha), 
we include two orbitals per unit cell labelled as m = A 
and m = B (the A atoms and the B atoms are each 
arranged on an hexagonal Bravais lattice). 

We begin by considering the diagonal element Haa 
of the transfer integral matrix H, equation ([6]), for the 
A site orbital. It may be determined by substituting the 
Bloch function ([3]) for m — A into the matrix element (|6|) , 
which results in a double sum over the positions of the 
unit cells in the crystal. Assuming that the dominant 
contribution arises from those terms involving a given 
orbital interacting with itself (i.e., in the same unit cell), 
the matrix element may be written as 



Haa 



1 N 

-Y 



(r-R Ai i)|^|^(r-R A)i )). (8) 



»=i 



This may be regarded as a summation over all unit cells of 
a parameter ea = {4>a (r — Ra,i) |W|<Aa ( r — Ra,<)) that 
takes the same value in every unit cell. Thus, the matrix 
element may be simply expressed as Haa ~ ^A- Simi- 
larly, the diagonal element Hbb for the B site orbital can 
be written as Hbb = e_B, while for intrinsic graphene ca 
is equal to es as the two sublattices are identical. The 
calculation of the diagonal elements of the overlap inte- 
gral matrix S, equation ©j proceeds in the same way 
as that of H, with the overlap of an orbital with itself 



equal to unity, ( 
Sbb = Saa = 1- 



(r - \4>j (r - Rj,i)) = 1. Thus, 



The off-diagonal element Hab of the transfer integral 
matrix H describes the possibility of hopping between 
orbitals on A and B sites. Substituting the Bloch func- 
tion ^ into the matrix element (J5]) results in a sum over 
all A sites and a sum over all B sites. We assume that 
the dominant contribution arises from hopping between 
adjacent sites. If we consider a given A site, say, then we 
take into account the possibility of hopping to its three 
nearest-neighbour B sites, labelled by index I = 1,2,3: 

1 N 3 
i=l 1 = 1 

x (cf> A (r - R Ati ) \H\4> B (r - R A i - Si)) , (9) 

where 5i are the positions of three nearest B atoms rel- 
ative to a given A atom, which may be written as S\ = 
(0,a/V3), 8 2 = (a/2,-a/2V3), <5 3 = (-a/2, -a/2-v/t) . 

The sum with respect to the three nearest-neighbour 
B sites is identical for every A site. A hopping parameter 
may be defined as 



7o 



(r-RA,i)|H|Mr-RA,i-*i)>> ( 10 ) 



which is positive. Then, the matrix element may be writ- 
ten as 



H 



AB 



-70/ (k); /(k) 



i=i 



(11) 



The other off-diagonal matrix element is given by Hba = 
H*ab ~ ~7o/*(k). The function / (k) describing 
nearest- neighbor hopping, equation (ITTj) . is given by 

/ (k) = e zk v a/V * + 2e- lk * a l 2s/l cos (k x a/2) , (12) 

where k = [k x , k y ) is the in-plane wave vector. The cal- 
culation of the off-diagonal elements of the overlap in- 
tegral matrix S is similar to those of H. A parameter 
so = (4>a (r — Ra,») \4>b (r — Rs,z)) is introduced to de- 
scribe the possibility of non-zero overlap between orbitals 
on adjacent sites, giving Sab = S B a — s of (k). 

Gathering the matrix elements, the transfer H m and 
overlap S m integral matrices of monolayer graphene may 
be written as 



H m = 



ca -7o/ (k) 
-7o/* (k) e B 

1 so/ (k) 
sof* (k) 1 



(13) 
(14) 



The corresponding energy may be determined [531 ] by 
solving the secular equation (JT]). For intrinsic graphene, 
i.e., €a — (b = 0, we have 



E± 



±7ol/(k)| 
lTSo|/(k) 



(15) 



The parameter values are listed by Saito et al [53| as 
70 = 3.033 eV and s = 0.129. 

The function /(k), equation (fl"2"|) is zero at the corners 
of the Brillouin zone, two of which are non-equivalent 
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(i.e., they are not connected by a reciprocal lattice vec- 
tor). For example, corners K + and K_ with wave vectors 
K± = ±(47r/3a, 0) are labelled in FigureQJb). Such po- 
sitions are called K points or valleys, and we will use a 
valley index £ = ±1 to distinguish points JQ. At these 
positions, the solutions (fT5|) are degenerate, marking a 
crossing point and zero band gap between the conduction 
and valence bands. The transfer matrix H m is approxi- 
mately equal to a Dirac-like Hamiltonian in the vicinity 
of the K point, describing massless chiral quasiparticlcs 
with a linear dispersion relation. These points are partic- 
ularly important because the Fermi level is located near 
them in pristine graphene. 



3. Bilayer graphene 

In the tight-binding description of bilayer graphene, we 
take into account 2p z orbitals on the four atomic sites in 
the unit cell, labelled as j = A1,B1,A2,B2. Then the 
transfer integral matrix of bilayer graphene [9l [iTl 154 l6lT- 
l63j is a 4 x 4 matrix given by 



fli 



cai -70/ (k) 74/ (k) -73/* (k) 

. "To/* ( k ) e Bi 7i 74/ (k) 

1,1 74/* (k) 71 tA2 -70 f (k) 

-73/ (k) 74/* (k) -70/* (k) e B2 



(16) 

where the tight-binding parameters are defined as 

70 = -{<t>Ai\U\<j>Bx) = -<&u|H|(M, ( 17 ) 

71 = {(t>A2\H\<f>Bx), (18) 

l3 = -(<f>Ai\H\cl> B2 ), (19) 

74 = {^Al\n\^A2) = (4>Bl\n\(t>B2). (20) 

Here, we use the notation of the Slonczewski-Weiss- 
McClure (SWM) model [IHS^I that describes bulk 
graphite. Note that definitions of the parameters used 
by authors can differ, particularly with respect to signs. 

The upper-left and lower-right square 2x2 blocks of 
describe intra-layer terms and are simple generalisations 
of the monolayer, equation (|13p . For bilayer graphene, 
however, we include parameters describing the on-site 
energies e A i, e B \, e A2 , e B2 on the four atomic sites, that 
are not equal in the most general case. As there are four 
sites, differences between them are described by three 
parameters [63l |: 

*ai = \{-U + 8 AB ) , (21) 

em = § (-U + 2A' - 5 AB ) , (22) 

e A 2 = \{U + 2A' + S AB ) , (23) 

em = \{U-5 AB ) , (24) 

where 



U 
A' 



| [(em 



em) 
- £.42) 



{(-A2 + em)] , 

(cai + em)} , 



Sab = \ [(cai + £A2) - (em + em)} 



(25) 
(26) 
(27) 



The three independent parameters are U to describe in- 
terlayer asymmetry between the two layers 0, [l(J 
[20l . E3-[z3, A' for an ener gy d ifference between dimer 
and non-dimer sites (5^ - [5a [67| . and 5ab for an energy 
difference between A and B sites on each layer [63j, l75j . 
These parameters are described in detail in sections III Fl 
and |ntO 

The upper-right and lower-left square 2x2 blocks of 
.ffb describe inter-layer coupling. Parameter 71 describes 
coupling between pairs of orbitals on the dimer sites Bl 
and A2: since this is a vertical coupling, the correspond- 
ing terms in H h (i.e., H A 2,bi = H B i,A2 = 7i) do not 
contain / (k) which describes in-plane hopping. Param- 
eter 73 describes interlayer coupling between non-dimer 
orbitals Al and B2, and 74 describes interlayer coupling 
between dimer and non-dimer orbitals Al and A2 or Bl 
and B2. Both 73 and 74 couplings are 'skew': they are 
not strictly vertical, but involve a component of in-plane 
hopping, and each atom on one layer (e.g., Al for 73) has 
three equidistant nearest-neighbours (e.g., B2 for 73) on 
the other layer. In fact, the in-plane component of this 
skew hopping is analogous to nearest-neighbour hopping 
within a single layer, as parameterised by 70- Hence, 
the skew interlayer hopping (e.g., H A i B2 — —73/* (k)) 
contains the factor / (k) describing in-plane hopping. 

It is possible to introduce an overlap integral matrix 
for bilayer graphene |63| 



Sb 



( 1 s /(k) 

s /*(k) 1 Sl 

si 1 s /(k) 

s /*(k) 1 



,(28) 



with a form that mirrors H^. Here, we only include 
two parameters: s = (0ai|^bi) = (4>A2\4>B2} describing 
non-orthogonality of intra-layer nearest-neighbours and 
Si = ((/^I^bi) describing non-orthogonality of orbitals 
on dimer sites Al and B2. In principle, it is possible 
to introduce additional parameters analogous to 73, 74, 
etc., but generally they will be small and irrelevant. In 
fact, it is common practice to neglect the overlap inte- 
gral matrix entirely, i.e., replace Sb with a unit matrix, 
because the influence of parameters sq and S\ describ- 
ing non-orthogonality of adjacent orbitals is small at low 
energy \E\ < 71. Then, the generalised eigenvalue equa- 
tion ([5]) reduces to an eigenvalue equation H^ip = Eip 
with Hamiltonian H^, equation (1161) . 

The energy differences U and S AB are usually at- 
tributed to extrinsic factors such as gates, substrates 
or doping. Thus, there are five independent parame- 
ters in the Hamiltonian (|16|) of intrinsic bilayer graphene, 
namely 70, 71, 73, 74 and A'. The band structure pre- 
dicted by the tight-binding model has been compared to 
observations from photoemission 1161 . Raman |76| and 
infrared spectroscopy [H, [56|, [78l48ll ]. Parameter values 
determined by fitting to experiments are listed in Table U 
for bulk graphite [67j|, for bilayer graphene by Raman 
[zl, [73] and infrared [H, [H, H3] spectroscopy, and for 
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TABLE I. Values (in eV) of the Slonczewski-Weiss-McClure (SWM) model parameters [&il467| determined experimentally. 
Numbers in parenthesis indicate estimated accuracy of the final digit (s). The energy difference between dimer and non-dimer 
sites in the bilayer is A' = A — 72 +75. Note that next-nearest layer parameters 72 and 75 are not present in bilayer graphene. 



Parameter Graphite [67] Bilayer [76] Bilayer [55] Bilayer [56] Bilayer [80] Trilayer [82] 

70 3.16(5) 2.9 3.(f - 3.16(3) 3.1®. 

71 0.39(1) 0.30 0.40(1) 0.404(10) 0.381(3) 0.39^ 

72 -0.020(2) - -0.028(4) 

73 0.315(15) 0.10 0.3.1 " 0.38(6) 0.315^ 

74 0.044(24) 0.12 0.15(4) - 0.14(3) 0.041(10) 

75 0.038(5) - 0.05(2) 
A -0.008(2) - 0.018(3) 0.018(2) 0.022(3) -0.03(2) 
/V 0.050(6) - 0.018(3) 0.018(2) 0.022(3) 0.046(10) 



a This parameter was not determined by the given experiment, the value quoted was taken from previous literature. 



Bernal-stacked trilayer graphene by observation of Lan- 
dau level crossings [82| . Note that there are seven param- 
eters in the Slonczewski-Weiss-McClure (SWM) model of 
graphite |64| - |67| because the next-nearest layer couplings 
72 and 75, absent in bilayer, are present in graphite (and 
trilayer graphene, too). Parameter A in the SWM model 
is related by A = A'+72— 75 to the parameter A' describ- 
ing the energy difference between dimer and non-dimer 
sites in bilayer graphene. 

The energy bands are plotted in figure [3] along the 
k x axis in reciprocal space intersecting the corners K_, 
K + and the centre T of the Brillouin zone [see fig- 
ure IHb)]. Plots were made using Hamiltonian i?b, equa- 
tion (1161) . with parameter values determined by infrared 
spectroscopy 70 = 3.16 eV, 71 = 0.381 eV, 73 = 0.38 eV, 
74 = 0.14eV, e B i = e A 2 = A' = 0.022eV, and e A i = 
(B2 = U = Sab = [80]. There are four bands because 
the model takes into account one 2p z orbital on each of 
the four atomic sites in the unit cell; a pair of conduc- 
tion bands and a pair of valence bands. Over most of 
the Brillouin zone, each pair is split by an en ergy of the 
order of the interlayer spacing 71 w 0.4 eV [83] ]. Near 
the K points, inset of figure G2 one conduction band and 
one valence band are split away from zero energy by an 
energy of the order of the interlayer coupling 71 , whereas 
two bands touch at zero energy [jj]. The 'split' bands are 
a bonding and anti-bonding pair arising from the strong 
coupling (by interlayer coupling 71 ) of the orbitals on the 
dimer Bl and A2 sites, whereas the 'low-energy' bands 
arise from hopping between the non-dimer Al and B2 
sites. In pristine bilayer graphene, the Fermi level lies at 
the point where the two low-energy bands touch (shown 
as zero energy in figure [3]) and, thus, this region is rele- 
vant for the study of electronic properties. It will be the 
focus of the following sections. 

At low energy, the shape of the band structure pre- 
dicted by the tight-binding model (see inset in figure[3]) is 
in good agreement with that calculated by density func- 
tional theory 0, [5?], H|[ and it is possible obtain values 
for the tight-binding parameters in this way, generally 




FIG. 3. Low-energy bands of bilayer graphene arising from 
2p z orbitals plotted along the k x axis in reciprocal space in- 
tersecting the corners K-, K+ and the centre T of the Bril- 
louin zone. The inset shows the bands in the vicinity of the 
K+ point. Plots were made using Hamiltonian Hb, equa- 
tion (|16p . with parameter values 70 = 3.16 eV, 71 = 0.381 eV, 
73 = 0.38 eV, 74 = 0.14 eV-em = e A 2 = A' = 0.022 eV, and 
€ A i = e B2 = U = 8 AB = [IJ. 



in line with the experimental ones listed in Table [I] The 
tight-binding model Hamiltonian Hh , equation (IT6"|) , used 
in conjuction with the parameters listed in Table [H is not 
accurate over the whole Brillouin zone because the fitting 
of tight-binding parameters is generally done in the vicin- 
ity of the corners of the Brillouin zone K + and K- (as 
the Fermi level lies near zero energy). For example, pa- 
rameter So in equation l|28p describing non-orthogonality 
of adjacent orbitals has been neglected here, but it con- 
tributes electron-hole asymmetry which is particularly 
prevalent near the T point at the centre of the Brillouin 
zone HIH. 
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4- Effective four-band Hamiltonian near the K points 

To describe the properties of electrons in the vicinity of 
the K points, a momentum p = frk — KK^ is introduced 
which is measured from the centre of the Kt point. Ex- 
panding in powers of p, the function /(k), equation (1T21) . 
is approximately given by /(k) w — V3a(£,p x — ip y )/2h 
which is valid close to the point, i.e., for pa/h <C 1, 
where p = |p| = (p 2 + Pp 1 ^ 2 - In monolayer graphene, 
the Hamiltonian matrix (|13|) is simplified by keeping only 
linear terms in momentum p as 



Hm = 



V7T €b 



(29) 



where 7r = t;p x + ip y , = £p x — ip y , and v = V3a~/o/2h 
is the band velocity. In the intrinsic case, 6a = £b = 0, 
the eigen energy becomes E = ±u|p|, which approxi- 
mates Eq. (|15[) . In bilayer graphene, similarly, Eq. (|16[) 
is reduced to 



H h = 







— 


V 3 TT 


VTT 




7i 


— V^TT 




7i 


£.42 




V3ir^ 


— «47T 


VTT 


£B2 



(30) 



where we introduced the effective velocities, v 3 = 
s/3a"fa/2H and V4 = s/3a"f4,/2h. 

At zero magnetic field Hamiltonian (|30|) yields four 
valley-degenerate bands -E'(p). A simple analytic solution 
may be obtained by neglecting the terms v^k, v^-k* pro- 
portional to 74, and by considering only interlayer asym- 
metry U in the on-site energies: £,41 = esi = —U/2 and 
e.A2 = (B2 = U/2. Then, there is electron-hole symme- 
try, i.e., energies may be written E = ±e Q (p), a = 1,2, 
H with 



0? 



U 2 
4 



P 



(-i) a VT, (3i) 



1/ 2 2 2\ 2 1 2 2 r 2 1 TT 2 . 2 21 

r = 3 (7i - v 3 p ) + v p [jx+lJ +v 3 p \ 
+ 2£jiv 3 v 2 p 3 cos 3tp , 

where tp is the polar angle of momentum p = (jp x ,p v ) — 
p (cos (p, sin ip). Energy £2 describes the higher-energy 
bands split from zero energy by the interlayer coupling 
71 between the orbitals on the dimer sites Bl, A2. 

Low-energy bands E — ±£1 are related to orbitals on 
the non-dimer sites Al, B2. In an intermediate energy 
range U, (v 3 /v) 2 ji < e\ < 71 it is possible to neglect 
the interlayer asymmetry U and terms proportional to 
73 (i.e., set U = v% = 0), and the low-energy bands may 
be approximated [9j as 



hi 



1 + 4wV/7i 2 - 1 



(32) 



which interpolates between an approximately linear dis- 
persion £1 sa vp at large momentum to a quadratic one 



E\ ~ p 2 /2m at small momentum, where the mass is 
rn = ji/2v 2 (see inset in figure [3]). This crossover oc- 
curs at p f=a 7i/2w. A convenient way to describe the 
bilayer at low energy and momentum p -C 71 /2v is to 
eliminate the components in the Hamiltonian (J30J) re- 
lated to the orbitals on dimer sites Bl, A2, resulting in 
an effective two-component Hamiltonian describing the 
orbitals on the non-dimer sites Al, B2, and, thus, the 
two bands that approach each other at zero energy. This 
is described in the next section, and the solutions of this 
Hamiltonian are shown to be massive chiral quasiparti- 
cles [i|,|9j], as opposed to massless chiral ones in monolayer 
graphene. 



C. Effective two-band Hamiltonian at low energy 

In this section we focus on the low-energy electronic 
band structure in the vicinity of the points K+ and A"_ 
at the corners of the first Brillouin zone, relevant for en- 
ergies near the Fermi level. A simple model may be ob- 
tained by eliminating orbitals related to the dimer sites, 
resulting in an effective Hamiltonian for the low-energy 
orbitals. First, we outline the procedure in general terms, 
because it may be applied to systems other than bilayer 
graphene such as ABC-stacked (rhombohcdral) graphene 
multilayers f&il . f85j | , before applying it specifically to bi- 
layer graphene. 



1. General procedure 

We consider the energy eigenvalue equation, and con- 
sider separate blocks in the Hamiltonian corresponding to 
low-energy 9 = {ip A i,^B2) T and dimer x = (^Ai^Bif 
components: 



hg u 




= E 



X 



(33) 



The second row of (|33|) allows the dimer components to 
be expressed in terms of the low-energy ones: 



x = (E-K 



(34) 



Substituting this into the first row of (|33|) gives an effec- 
tive eigenvalue equation written solely for the low-energy 
components: 



-v 



6 = E6. 



9 + u(E - , 

[h e - uh' 1 ^] 6 « ESQ , 

where 5=1 + uh^ 2 ^ . The second equation is accurate 
up to linear terms in E. Finally, we perform a transfor- 
mation $ = S 1/2 9: 



[h - u/i" V] S- 1/2 $ « ES 1/2 1> , 
S~ 1/2 [h e - uh- l v)] S- 1/2 1> sa . 



(35) 
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This transformation ensures that normalisation of $ is 
consistent with that of the original states: 

$t$ = 0ts0 = f (1 + uh- 2 v)) 9 , 

~e f e + x f x, 

where we used equation p4[) for small E: x ~ — h^u^O. 
Thus, the effective Hamiltonian for low-energy compo- 
nents is given by equation ()35|) : 



ff (eff) w 5 -l/2 _ u/l -l u t] 5 -l/2 f (36) 

S=l + uhZ 2 v). (37) 



2. Bilayer graphene 



The Hamiltonian (1301) is written in basis 
Al,Bl,A2,B2. If, instead, it is written in 
the basis of low-energy and dimer components 
(0,x) = A1,B2,A2,B1, equation ([33]), then 



e A i v 3 tt 

V 3 TT^ e B 2 



— t>47P 



VTT —V^TT 



£A2 71 

7i e-Bi / 
/ — v^ir 

\ VTT 



Using the procedure described in the previous section, 
equations (|36I37[) . it is possible to obtain an effective 
Hamiltonian _ff( cff ) = H 2 for components (i/jai, 1 pB2)- An 
expansion is performed by assuming that the intralayer 
hopping 70 and the interlayer coupling 71 are larger than 
other energies: 70,71 > \E\, vp, \<y 3 \, \<y 4 \, \U\, \A'\, \5 AB \- 
Then, keeping only terms that are linear in the small 
parameters I73I, I74I, \U\, |A'|, \5ab\ and quadratic in mo- 
mentum, the effective Hamiltonian [H, |63| is 



'-'3 



H2 = ho + h w + hi + Iia + hu + \iab , 
, l/0 (^t) 2 \ 

/ 7T \ _ v 3 a ( (tt+) : 
?rt J 4V3h \tt 2 

2iw 4 I tt^it 

71 \ 7T7rt 
AV I TT^TT 



(38) 



ITTT^ 



fin = - 



hAB — 




2v 2 I t^tt 

7M 



where n = £p x + ip y , ir^ — £p x — ip y . In the following 
sections, we discuss the terms in Hi- The first term ho 



describes massive chiral electrons, section Hi Dl It gener- 
ally dominates at low energy \E\ <§; 71, so that the other 
terms in H2 may be considered as perturbations of it. 
The second term h w , section HTeI introduces a triangu- 
lar distortion of the Fermi circle around each K point 
known as 'trigonal warping'. Terms hu and hAB, with 
±1 on the diagonal, produce a band gap between the 
conduction and valence bands, section III G[ whereas 
and ft.A introduce electron-hole asymmetry into the band 
structure, section lITFl 

The Hamiltonian (|38j) is written in the vicinity of a 
valley with index £ = ±1 distinguishing between K + and 
A'_ . In order to briefly discuss the effect of symmetry op- 
erations on it, we introduce Pauli spin matrices cr x , a y , a z 
in the A1/B2 sublattice space and H x , H y , H z in the val- 
ley space. Then, the first term in the Hamiltonian may be 
written as h Q = -(l/2m)[a x (pl~pl) + 2U z a y p x p y ]. The 
operation of spatial inversion i is represented by H x a x 
because it swaps both valleys and lattice sites, time inver- 
sion is given by complex conjugation and n z , as it swaps 
valleys, too. Hamiltonian (|35|) satisfies time-inversion 
symmetry at zero magnetic field. The intrinsic terms 
ho, h w , /14, and /ia satisfy spatial-inversion symmetry 
because the bilayer crystal structure is spatial-inversion 
symmetric, but terms hjj and hAB, with ±1 on the di- 
agonal, are imposed by external fields and they violate 
spatial-inversion symmetry, producing a band gap be- 
tween the conduction and valence bands. 



D. Interlayer coupling 71: massive chiral electrons 



The Hamiltonian ho in equation (|38|) resembles the 
Dirac-like Hamiltonian of monolayer graphene, but 
with a quadratic-in-momentum term on the off-diagonal 
rather than linear. For example, the term it 2 /2m ac- 
counts for an effective hopping between the non-dimer 
sites Al, B2 via the dimer sites Bl, A2 consisting of a 
hop from Al to Bl (contributing a factor vtt), followed by 
a transition between Bl, A2 dimer sites (giving a 'mass' 
~ 71), and a hop from A2 to B2 (a second factor of vtt). 
The solutions are massive chiral electrons @, with 
parabolic dispersion E = ±p 2 /2m, m = ji/2v 2 . The 
density of states is m/(2irh 2 ) per spin and per valley, and 
the Fermi velocity vf — Pf/tti is momentum dependent, 
unlike the Fermi velocity v of monolayer graphene. 

The corresponding wave function is given by 



-M 1 

V2 \ Tt 



ip.r /h 



(39) 



The wave function components describe the electronic 
amplitudes on the Al and B2 sites, and it can be useful 
to introduce the concept of a pseudospin degree of free- 
dom [1, Q that is related to these amplitudes. If all the 
electronic density were located on the Al sites, then the 
pseudospin part of the wave function |f) = (1,0) could be 
viewed as a pseudospin 'up' state, pointing upwards out 
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of the graphene plane. Likewise, a state \l) = (0, 1) with 
density solely on the B2 sites could be viewed as a pseu- 
dospin 'down' state. However, density is usually shared 
equally between the two layers, so that the pseudospin 



is a linear combination of up and down, \~\) =p 



II), 



equation |39|) . and it lies in the graphene plane. 

The Hamiltonian may also be written as ho — 
{p 2 /2m) <r.n 2 where the pseudospin vector is er = 
(a x , <j y ,<j z ), and n.2 = — (cos 2ip, £ sin 2(p, 0) is a unit vec- 
tor. This illustrates the chiral nature of the electrons 
0,0: the chiral operator er.n 2 projects the pseudospin 
onto the direction of quantization n.2 , which is fixed to lie 
in the graphene plane, but turns twice as quickly as the 
momentum p. For these chiral quasiparticles, adiabatic 
propagation along a closed orbit produces a Berry's phase 
[86| change of 2n [1, Q of the wave function, in contrast 
to Berry phase tt in monolayer graphene. 

Note that the chiral Hamiltonian h may be viewed as a 
generalisation of the Dirac-like Hamiltonian of monolayer 
graphene and the second (after the monolayer) in a family 
of chiral Hamiltonians Hj, J = 1,2,..., corresponding to 
Berry's phase Jn which appear at low energy in ABC- 
stacked (rhombohedral) multilayer graphene (914151 HH, 



Hj 



gj 



(7T+)* 
7T J 



(40) 



where g\ = v for monolayer, g 2 = — l/2m for bilayer, and 
.93 = v 111 f° r trilayer graphene. Since the pseudospin is 
related to the wavefunction amplitude on sites that are 
located on different layers, pseudospin may be viewed as 
a 'which layer' degree of freedom fLU l87j|. 



Interlayer coupling 73: trigonal warping and the 
Lifshitz transition 



E. 



The Hamiltonian ho in equation (|38|) yields a 
quadratic, isotropic dispersion relation E — ±p 2 /2m 
with circular iso-energetic lines, i.e., there is a circular 
Fermi line around each K point. This is valid near the 
K point, pa/h <C 1, whereas, at high energy, and mo- 
mentum p far from the K point, there is a triangular 
perturbation of the circular iso-energetic lines known as 
trigonal warping, as in monolayer graphene and graphite. 
It occurs because the band structure follows the symme- 
try of the crystal lattice as described by the full momen- 
tum dependence of the function /(k), equation (fT?l) [88l . 

In bilayer graphene @, as in bulk graphite |89l - l92| . a 
second source of trigonal warping arises from the skew in- 
terlayer coupling 73 between non-dimer A\ and B2 sites. 
The influence of 73 on the band structure is described 
by equation (|3ip . In the two-band Hamiltonian, it is 
described by h w in equation (|38p . the second term of 
which arises from a quadratic term in the expansion of 
/(k) a -V5o(CPx - ip y )/2fr + a 2 {& x + i Py ) 2 /8h 2 . This 
second term has the same momentum dependence as ho, 



and, thus, it actually only gives a small additional con- 
tribution to the mass m. The first term in h w causes 
trigonal warping of the iso-energetic lines in directions 
ip = <po, where ipo — 0, |7T, ^n at K + , <p = ^n,n, |7r at 
K . 

To analyse the influence of h w at low energy, we con- 
sider just ho and the first term in h w , and the resulting 
energy E — ±e± is given by 



(v 3 p) 



£v 3 p 3 



cos {3(f) + ( 

1 2m 



(41) 



in agreement with equation pip for [7 = 0, wp/71 <C 
1, and »3/o < 1. As it is linear in momentum, the 
influence of h w and the resulting triangular distortion 
of iso-energetic lines tend to increase as the energy and 
momentum are decreased until a Lifshitz transition [93[ 
occurs at energy 



71 / W3 \ ' 
4 U 



1 meV. 



(42) 



For energies \E\ < sl, iso-energetic lines are broken into 
four separate 'pockets' consisting of one central pocket 
and three 'leg' pockets, the latter centred at momentum 
p I1V3/V 2 and angle tpo, as shown in Figure 2J The 
central pocket is approximately circular for \E\ <C £l 
with area A c ~ ire 2 / (hv 3 ) 2 , while each leg pocket is ap- 
proximately elliptical with area Ai ~ A c /3- Note that 
Berry phase 2ir is conserved through the Lifshitz transi- 
tion; the three leg pockets each have Berry phase tt while 
the central pocket has — tt fl2l. l94l|. 



F. Interlayer coupling 74 and on-site parameter A': 
electron-hole asymmetry 

Skew interlayer coupling 74 between a non-dimer and 
a dimer site, i.e., between Al and A2 sites or between 
Bl and B2 sites, is described by /14 in equation (|38l) . 
where the effective velocity is V4 = v / 3i74/2?i. This term 
produces electron-hole asymmetry in the band struc- 
ture, as illustrated by considering the energy eigenvalues 
E = ±(p 2 /2m)(l ± 2v4/v) of the Hamiltonian ho + h,4- 
The energy difference A' between dimer and non-dimer 
sites, e^i = £B2 = 0, €bi — eA2 = A', equation (|26j) . 
also introduces electron-hole asymmetry into the band 
structure: the low-energy bands described by ho + hj\ 
are given by E = ±p 2 /2m(l ± A'/^i). 



G. Asymmetry between on-site energies: band 
gaps 

1. Interlayer asymmetry 

Interlayer asymmetry U , equation (|25p . describes a dif- 
ference in the on-site energies of the orbitals on the two 
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E/e L 



Px'Pl 



FIG. 4. (a) Trigonal warping of the equi-energy lines in the 
vicinity of each K point, and the Lifshitz transition in bilayer 
graphene. The energy is in units of sl- (b) Corresponding 
three-dimensional plot of the low-energy dispersion. 



2. Intralayer asymmetry between A and B sites 

The energy difference 5ab between A and B sites may 
be described by the Hamiltonian (|30|) with cai — ~ est = 
£^42 = — ( B2 — Sab/^ and v 3 = v 4 = 0, yielding bands 
E = ±e a : 



S AB 



^/l + 4^V/7i + (-1)" 



(45) 



Thus, Sab creates a band gap, but there is no Mexican 
hat structure. 



H. Next-nearest neighbour hopping 



The terms described in Hamiltonians (| 16I30I38[) do not 
represent an exhaustive list of all possibilities. Addi- 
tional coupling parameters may be taken into account. 
For exam ple, next-nearest neighbour hopping within each 
layer |95T - [9a ] results in a term (3 — |/(k)| 2 )7„ appear- 
ing on every diagonal element of the Hamiltonian Q16[). 
where 7„ is the coupling parameter between next-nearest 
A (or B sites) on each layer. Ignoring the constant-in- 
momentum part S"f n produces an additional term in the 
two-component Hamiltonian 



h„ = 



-y n v 2 p 2 (10 



7o 2 



1 



resulting in energies E = ±p 2 /2m(l =p 7n7i/7o)- 
Thus, next-nearest neighbour hopping represents another 
source of electron-hole asymmetry, after h 4l h&, and Sh- 



layers cai — est = — eA2 = — e_B2 = —U/2. Its influence 
on the bands E — ±e Q (p) is described by equation (|3"Tj) 
with U3 = 0: 



7? 



U 



7i 
4 



■U 2 ], 
(43) 



The low-energy bands, a = 1, display a distinctive 'Mex- 
ican hat' shape with a band gap U g between the conduc- 
tion and valence bands which occurs at momentum p g 
from the centre of the valley: 



U„ = 



\U\-n 



VlT+U 2 



\U\ /2t[+L72 
7l 2 + U 2 



(44) 



For small values of the interlayer asymmetry U, the band 
gap is equal to the asymmetry U g — \U\, but for large 
asymmetry values \U\ ^> 71 the band gap saturates 
U g — > 71. It is possible to induce interlayer asymme- 



try in bilayer graphene through doping [16| or the use of 
external gal 
section IIII1 



graphene th 

external gates [13, EH, H3 ■ This is described in detail in 



I. Spin-orbit coupling 

For monolayer graphene, Kane and Mele (9^| employed 
a symmetry analysis to show that there are two distinct 
types of spin-orbit coupling at the corners K + and K_ of 
the Brillouin zone. These two types of spin-orbit coupling 
exist in bilayer graphene, too. In both monolayers and 
bilayers, the magnitude of spin-orbit coupling - although 
the subject of theoretical debate - is generally considered 
to be very small, with estimates roughly in the range of 
1 to lOO/xeV [H4103. 



At the corner of the Brillouin zone in bilayer 
graphene, the contribution of spin-orbit coupling to the 
two-component low-energy Hamiltonian (I38[) may be 
written as 

hso = \so£,OzS z , (46) 

h.R = \r. (£<T X Sy + (JyS x ) , (47) 

where cr,, i = x,y, z are Pauli spin matrices in the Al / B2 
sublattice space, and Sj, j = x,y,z are Pauli spin ma- 
trices in the spin space. The first term hso is intrin- 
sic to graphene, i.e. it is a full invariant of the sys- 
tem. Both intra- and inter-layer contributions to hso 
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have been discussed [I05l4l07i Il09| with the dominant 
contribution to its magnitude Ago attribute d to skew 
interlayer coupling between n and a orbitals |l06l . Il07j | 
or to the presence of unoccupied d orbitals within each 
graphene layer |l 09j j . Taken with the quadratic term ho 
in the Hamiltonian ([55]) , h so produces a gap of magni- 
tude 2A50 in the spectrum of bilayer graphene, but the 
two low-energy bands remain spin and valley degenerate 
(as in a monolayer): E — ±yj A| Q + v 4 p 4 /jf. However, 
there are gap less edge excitations and, like monolayer 
graphene [9Pj , bilayer graphene in the presence of intrin- 
sic spin-orbit coupling is a t opologica l insulator with a 
finite spin Hall conductivity [llOL llllj . 

The second type of spin-orbit coupling is the Bychkov- 
Rashba term hn, equation (|47[) . which is permitted only if 
mirror reflection symmetry with respect to the graphene 
plane is broken, by the pres e nce of an electric field or 
a substrate, say pl- flOli fTo l [TToMTTej . Taken with the 
quadratic term ho in the Hamiltonian (|38j) . hji does not 
produce a gap, but, as in the monolayer, spin-splitting 
of magnitude 2A^ between the bands. That is, there are 
four valley-degenerate bands at low energy, 



(48) 



Generally speaking, there is a rich interplay between 
tuneable interlayer asymmetry U and the influence of the 
intrinsic and the Bychkov-Rashba spin-or bit coupling in 
bilayer graphene $M EE El EE EE ■ For example, 
the presence of interlayer asymmetry U breaks inversion 
symmetry and allows for spin-split levels in t he pr esence 
of intrinsic spin-orbit coupling only (Ai? = 0) [105| . while 
the combination of finite U and very large Rashba cou- 
pling has been predicted t o lea d to a topological insulator 
state even with Xso = EE. 




J. The integer quantum Hall effect 

When a two-dimensional electron gas is placed in a 
perpendicular magnetic field, electrons follow cyclotron 
orbit s and their energies are quantised as Landau levels 
[l!7j . At a high enough magnetic field strength, the dis- 
crete nature of the Landau level spectrum is manifest as 
the integer quantum Hall effect jll8rtl20j |. whereby the 
Hall conductivity assumes values that are integer multi- 
ples of the quantum of conductivity e 2 /h. 

The Landau level spe ctrum of monolayer graphene was 
calculated by McClure [121] nearly sixty years ago, and 
there hav e been a number of related theoretical studies 
[E IE HE Il22rtl24j considering the consequences of chi- 
rality in graphene. The experimental observation of the 
integer quantum Hall effect in monolayer graphene [E Q 
found an unusual sequencing of the quantised plateaus 
of Hall conductivity, confirming the chiral nature of the 
electrons and prompting an explosion of interest in the 



field Q. In bilayer graphene, the observation of the in- 
teger quantum Hall effect [8| and the calculated Landau 
level spectrum @ uncovered additional features related 
to the chiral nature of the electrons. 



1. The Landau level spectrum of bilayer graphene 

We consider the Landau level spectrum of the two- 
component chiral Hamiltonian ho, equation (|38|) . The 
magnetic field is accounted for by the operator p = 
(Px,Py) = —ihV+eA where A is the vector potential and 
the charge of the electron is — e. For a magnetic field per- 
pendicular to the bilayer, B = (0, 0, —B) where B = |B|, 
the vector potential may be written in the Landau gauge 
A = (0, — Bx, 0), which preserves translational invari- 
ance in the y direction. Then, tt = —iht;d x + hd y — ieBx 
and 7r ' = —iH£d x — hd y + ieBx, and eigenstates are com- 
prised of functions that are harmonic oscillator states 
in the x direction and plane waves in the y direction 

PIEE, 



4>i (x,y) = AiHel — 

As 



x exp 



Py^l 

h 

X 

As " 



Py^B 



+ i 



PyV 



,(49) 



where the magnetic length is A^ = \J h/eB, Tit are 
Hermite polynomials of order I for integer I > 0, and 
Ai = 1/ y^Zly^ is a normalisation constant. 

The operators tt and ir^ appearing in the Hamilto- 
nian (|38[) act as raising and lowering operators for the 
harmonic oscillator states (14T)|) . At the first valley, K + , 



y/2ih 



Vi + 1 <f>e+i 



(50) 
(51) 



and TT(f>o — 0. Then, it is possible to show that the Lan- 
dau level spectrum of the Hamiltonian (f3"51) consists of a 
series of electron and hole levels with energies and wave 
functions [i| given by 



K 



+ ■ 



E i>± = ±hu c ^e(£-l) , £ > 2, (52) 



where uj c — eB/m and ± refer to the electron and hole 
states, respectively. For high values of the index, £ ^> 1, 
the levels are approximately equidistant with spacing hio c 
proportional to the magnetic field strength B. However, 
this spectrum, equation (|52l) , is only valid for sufficiently 
small level index and magnetic field £Huj c <C 71 because 
the effective Hamiltonian (|38p is only applicable at low 
energy. 
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In addition to the field-dependent levels, there are two 
levels fixed at zero energy E\ = Eq = with cigcnfunc- 
tions: 



(a) CT xy (4eVh) 



K+ 



^ = 











,° , (54) 



They may be viewed as arising from the square of the 
lowering operator in the Hamiltonian (|38f which acts on 
both the oscillator ground state and the first excited state 
to give zero energy 7r 2 </>o = 7r 2 0i = 0. The eigenfunctions 
rpo and ipi have a finite amplitude on the Al sublattice, 
on the bottom layer, but zero amplitude on the B2 sub- 
lattice. 

At the second valley, K- , the roles of operators ir and 
7r^ are reversed: 



K- 



TT(pi 



7T t 0£ = 



y/2ih 



V2ih 



(55) 

Vlfa-i, (56) 



and 7r^0o = 0. The Landau level spectrum at K_ is 
degenerate with that at K + , i.e., Et t ± = ±huj c ^£(£ — 1) 
for £ > 2 and Ei = E = 0, but the roles of the Al and 
B2 sublattices are reversed: 



K. 



1 I 4>l-2 




£>2, (57) 



(58) 



The valley structure and electronic spin each contribute a 
twofold degeneracy to the Landau level spectrum. Thus, 
each level in bilayer level graphene is fourfold degenerate, 
except for the zero energy levels which have eightfold 
degeneracy due to valley, spin and the orbital degeneracy 
of 4> , ip%. 



2. Three types of integer quantum Hall effect 

The form of the Landau level spectrum is manifest in 
a measurement of the integer quantum Hall effect. Here, 
we will compare the density dependence of the Hall con- 
ductivity a xy (n) for bilayer graphene with that of a con- 
ventional semiconductor and of monolayer graphene. 

The Landau level spectrum of a conventional two- 
dimensional semiconductor is E# = haj c (£ + 1 /2), £ > , 
where lo c = eB/ra is the cyclotron frequency |l!9l . Il2(jj . 
As density is changed, there is a step in a xy whenever a 
Landau level is crossed, and the separation of steps on 
the density axis is equal to the maximum carrier den- 
sity per Landau level, gB/tpo, where ipo — h/e is the 
flux quantum and g is a degeneracy factor. Each plateau 
of the Hall conductivity a xy occurs at a quantised value 
of Nge 2 /h where N is an integer labelling the plateau 
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FIG. 5. Schematic of the dependence of the Hall conductiv- 
ity cr X y on carrier density n for (a) monolayer graphene and 
(b) bilayer graphene, where ipo — h/e is the flux quantum 
and B is the magnetic field strength. The dashed line in (b) 
shows the behaviour for a conventional semiconductor or bi- 
layer graphene in the presence of large interlayer asymmetry 
U (section III J 3|) with fourfold level degeneracy due to spin 
and valley degrees of freedom. 



and g is an integer describing the level degeneracy; steps 
between adjacent plateaus have height ge 2 /h. 

T he Land au level spectrum of monolayer graphene 
[98l Il2ll4l24 1 consists of an electron and a hole series of 
levels, Ei.± = ±^/2£hv / Xb for i > 1, with an additional 
level at zero energy Eq = 0. All of the levels are fourfold 
degenerate, due to spin and valley degrees of freedom. 
The corresponding Hall conductivity is shown schemati- 
cally in figure [5ja). There are steps of height 4e 2 /h be- 
tween each plateau, as expected by consideration of the 
conventional case, but the plateaus occur at half-integer 
values of 4e 2 //i instead of integer ones, as observed ex- 
perimentally @, Q . This is due to the existence of the 
fourfold-degenerate level Eq = at zero energy, which 
contributes to a step of height 4e 2 /h at zero density. 

For bilayer graphene, plateaus in the Hall conductiv- 
ity cF xy {n), Fig. [5lb), occur at integer multiples of 4e 2 /h. 
This is similar to a conventional semiconductor with level 
degeneracy g — 4 arising from the spin and valley de- 
grees of freedom. Deviation from the conventional case 
occurs at low density. In the bilayer there is a step in 
a xy of height 8e 2 /h across zero density, accompanied by 
a plateau separation of 8B/ipo in density [1,0], arising 
from the eightfold degeneracy of the zero-energy Lan- 
dau levels. This is shown as the solid line in Fig. [5jb), 
whereas, for a conventional semiconductor, there no step 
across zero density (the dashed line). 

Thus, the chirality of charge carriers in monolayer and 
bilayer graphene give rise to four- and eight-fold degen- 
erate Landau levels at zero energy and to steps of height 
of four and eight times the conductance quantum e 2 /h in 
the Hall conductivity at zero density [2j, y, |8( . Here, we 
have assumed that the degeneracy of the Landau levels 
is preserved, i.e., any splitting of the levels is negligible 
as compared to temperature and level broadening in an 
experiment. The role of electronic interactions in bilayer 
graphene is described in section IVIII1 while we discuss 
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the influence of interlayer asymmetry on the Landau level 
spectrum and integer quantum Hall effect in the next sec- 
tion. 



3. The role of interlayer asymmetry 



The Landau level states, equations (|53|57|) . have dif- 
ferent amplitudes on the lower (^41 sublattice) and up- 
per (B2 sublattice) layers, with the role of the sublattice 
sites swapped at the two valleys. Thus, interlayer asym- 
metry U as described by the effective Hamiltonian hjj, 
equation (|3"8)l , leads to a weak splitting of the valley de- 
generacy of the levels [§] : 



t;Uhuj c 

27i 



(59) 



Such splitting is prominent for the zero-energy states Q , 
equations (|54l58p . because they only have non-zero am- 
plitude on one of the layers, depending on the valley: 



E 



£JJtvjj c 
7i 



(60) 
(61) 



When the asymmetry is large enough, then the splitting 
U of the zero energy levels from each valley results in 
a sequence of quantum Hall plateaus at all integer val- 
ues of 4e 2 //i including a plateau at zero density [17[, as 
observed experimentally [20] - This behaviour is shown 
schematically as the dashed line in figure [S^b). The Lan- 
dau level spectrum in the presenc e of large interlayer 
asymmetry U has been calcu lated 125H128I ]. including 
an analysis of level crossings (l28| and a self-consistent 
calculati on of the spectrum in the presence of an exter- 
nal gate |l26l . Il27l | , generalising the zero-field procedure 
outlined in the next section. 



III. TUNEABLE BAND GAP 
A. Experiments 

A tuneable band gap in bilayer graphene was first ob- 
served with angle-resolved photoemission of epitaxial bi- 
layer graphene on silicon carbide 16], and the ability to 
control the gap was demonstrated by doping with potas- 
sium. Since then, the majority of experiments probing 
the band gap have used single or dual gate devices based 
on exfoliated bilayer graphene flakes [HI, [2(3] ■ The band 
gap has now been observed in a number of different ex- 
periments including photoemissio n [161. ma gneto t rans - 
port [13], infrared spectr osco py 
electronic co mpres sibility |l3lL 
spectroscopy [1331 ] , and transport 

The gap obs erved in o ptics flF 
up to 250meV[m,[IIiI, 
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FIG. 6. Schematic of bilayer graphene in the presence of ex- 
ternal gates located at x — —Lb and x — +Lt, with po- 
tentials Vt and Vt, which are separated from the bilayer by 
media of dielectric constants £{, and Et, respectively. The bi- 
layer is modelled as two parallel conducting plates positioned 
at x — —co/2 and +co/2, separated by a region of dielectric 
constant e r . The layers have charge densities <7i = —em and 
o~ 2 = —eri2 corresponding to layer number densities ni and 
7J2 • Charge densities Obo and o~to (dashed lines) arise from the 
presence of additional charged impurities. 



(as the gap should saturate at the value of the interlayer 
coupling 71L Tra nsport m easurements show insulating 
behaviour |19l. [3lL [l34T - |l39l ] . but, generally, not the huge 
suppression of conductivity expected for a gap of this 
mag nitude, and this has been attribute d to edge states 
11401 , the presence of diso rder (l4lMl43l or disorder and 



chiral charge carriers [144|. Broadly speaking, transport 
seems to occur through different mechanisms in diff eren t 
temperature regimes with thermal activation [3ll . Il34| - 
136l Il38l Il39j ] at high temper ature (ab o ve, r oughly, 2 to 
50 K) and variable- rang e 

MEMS 

1393 or nearest- 
neighbour hopping [I35I Il38l ] at low temperature. 



B. Hartree model of screening 

External gates are generally used to control the den- 
sity of electrons n on a graphene device [l[, but, in bi- 
layer graphene, external gates will also place the sepa- 
rate layers of the bilayer at different potential energies 
resulting in interlayer asymmetry U — £2 — £1 (where 
ei = e.4i = £si = -17/2 and £ 2 = £A2 = £52 = U/2). 
Thus, changing the applied gate voltage(s) will tend to 
tune both the density n and the interlayer asymmetry 
U, and, ultimately, the band gap U g . The dependence of 
the band gap on the density U g (n) relies on screening by 
electrons on the bilayer. In the follo wing, we d escribe a 
simple model [H M, H3, M, M, El EH that has 
been developed to take into account screening using the 
tight-binding model and Hartree theory. 



1. Electrostatics: asymmetry parameter, layer densities and 
external gates 

We use the SI system of units, and the electronic charge 
is — e where e > 0. The bilayer graphene device is mod- 
elled as two parallel conducting plates that are very nar- 
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row in the x-direction, continuous and infinite in the y- 
z plane, positioned at x — — c /2 and +cq/2, figure |U 
Here, Cq is the interlayer spacing and we denote the di- 
electric constant of the interlayer space as e r (it doesn't 
include the screening by 7r-band electrons that we are ex- 
plicitly modelling). Layer number densities are n\ and 
ri2, with corresponding charge densities a± = —en\ and 
02 = -en 2 . 

We take into account the presence of a back gate and a 
top gate, infinite in the y-z plane, located at x = —Lb and 
x = +L t , with potentials V& and V t , which are separated 
from the bilayer by media of dielectric constants £\, and 
e t , respectively. It is possible to describe the presence 
of additional charge near the bilayer - due to impurities, 
say - by introducing density nbo on the back-gate side 
and rito on the top-gate side. They correspond to charge 
densities <Jbo = en-bo and ato = erito, assuming that nbo 
and nto are positive for positive charge. 

Using Gauss's Law, it is possible to relate the external 
gate potentials and impurities concentrations t o the laye r 
densities and the interlayer asymmetry [Tt], Il45l 1147] : 

£Q£bVb EQStVt /„r,\ 

n = m + n 2 = — 1 h n b0 + n m , (62) 



eL h 



eL t 



e 2 c 



U=-^^eVt + — (n 2 -n t o) j ( 63 ) 

where the field within the bilayer interlayer space is ap- 
proximately equal to U/(eco). Equation (|62j) expresses 
the total electron density n = n\ + n 2 in terms of the ex- 
ternal potentials, generalising the relation for monolayer 
graphene [l[ . Note that the background densities nbo and 
nto shift the effective values of the gate potentials 14 and 
Vf. The second equation, for U, may be rewritten as 



U = U cxt + A 7l 
ec / £b 



(n 2 - ni) 



u '" = zrAf b v > 



St 



V, 



(64) 

A 7l (n&0 ~ nto) , (65) 



where the characteristic density scale nj_ and the dimen- 
sionless screening parameter A are 



7? 



irh 2 v 2 



A 



c e 2 7i 



(66) 



Equation (|64l) relates the asymmetry parameter U to a 
sum of its value, £/ cx t, if screening were negligible plus 
a term accounting for screening. Parameter values 71 = 
0.39eV and v = 1.0 x 10 6 ms" 1 give n ± = 1.1 x 10 13 cm" 2 . 
With interlayer spacing cq = 3. 35 A and an estimate for 
the dielectric constant of e r ~ 1, then A ~ 1, showing 
that screening is relevant. 



2. Calculation of individual layer densities 

Through electrostatics, the asymmetry parameter U 
is related to layer densities n\ and n 2 , equation 
The densities m and rii also depend on U because of 



the band structure of bilayer grap hene. Analytical cal- 
culations are possible 1?], 3 Il45j | if only the dominant 
inter-layer coupling 71 is taken into account in the four- 
band Hamiltonian (|3"0)) . Here we will use the two-band 
model (|38|) with an explicit ultraviolet cutoff [18| when 
integrating over the whole Brillouin zone. The simplified 
two-component Hamiltonian is 



Ho 



v i ( ° (^y 



71 








(67) 



Solutions to the energy eigenvalue equation Hii\) = Eip 
are given by 



E = ±a 



u 2 


w 4 n 4 




if 



E-U/2 



2E 



(68) 
(69) 



7i(B-(7/2)' 



Layer densities are determined by integration over the 
circular Fermi surface 



"-1(2) = ^2 / 1^1(52) (p) I pdp, (70) 

where a factor of four takes account of spin and valley 
degeneracy. 

For simplicity, we assume the Fermi energy lies within 
the conduction band. Using the solution (l69l) . the contri- 
bution of the part i ally- filled conduction band to the layer 
densities [17j, Il45l . Il48| is given by 



,ch 



(71) 



pdp 



^ rPF jj rPF 

rf? J Pdp T 2^ J 7^ 2 /4 + vV/7i 2 



n n±U ^ I 2|n|7i + ^j^ , ( 2 «7i 



2 4 7 i \n ± \U\ 



n±U 



where the minus (plus) sign is for the first (second) layer, 
Pf is the Fermi momentum, and the total density is n — 
p'p/irh 2 measured with respect to the charge neutrality 
point, i.e., we assume that the point of zero density is 
realised when the Fermi level lies at the crossing point of 
the conduction and valence bands. 

In addition, we take into account the contribution of 
the filled valence band n^ 2 ^ to the individual layer den- 
sities. Note that, as the asymmetry parameter U varies, 
the filled valence band does not contribute to any change 
in the total density n, but it does contribute to the dif- 
ference n\ —ri2. This may be obtained by integrating 
with respect to momentum as in equation (|71l) . but in- 
troducing an ultraviolet cutoff p max = equivalent 
to n max = n± [18J. Then, the contribution of the filled 
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valence band [U El HI EH, EM EH is given by 



n l(2) ~ =t - 



(72) 



Combining the contributions of equations (JTTJ) and (|72p . 

cb I „vb 
1(2) + n l(2)' 



the individual layer density, ni(2) = n u2) + n i(2) ' ' s gi ven 




Note that some calculations 148] include only the contri- 
bution (I7T1) of the partially-filled conduction band, others 
18] include only the filled valence band (f72|) . 




3. Self- consistent screening 

Substituting the layer density (I75|) into the expres- 
sion JM]) describing screening gives an expression [13, 
l60l . Il45j ] for the density-dependence of the asymmetry 
parameter U : 



U(n) 



where U cx t is the asymmetry in the absence of screen- 
ing (|d5|) . The extent of screening is described by the 
logarithmic term with argument depending on n and 
U , and a prefactor proportional to the screening pa- 
rameter A ~ 1 (as discussed earlier). A common ex- 
perimental setup, especially for exfoliated graphene on 
a silicon substrate [lj-d, Q , includes a single back gate. 
Figure [7] shows the density-dependence of the band gap 
U g {n) plotted as the back-gate voltage 14 is varied for a 
fixed top-gate voltage V t . In this case, the influence of 
the top-gate voltage V t may be absorbed into an effec- 
tive offset-density n = 2[eoe t Vt/(eL t ) + n m ] 17] giving 
U cx t = A7i(n — no)/nj_ in equation [TU Figure [5] shows 
the dependence of the difference in layer densities n\ — «2 
for the case no = including both the contribution of 
the partially filled bands as measured with respect to 
the charge neutrality point (fTTj) (dashed line) and the 
contribution of the full valence band (|72|) (dotted line). 
The sum of both terms (solid line) shows that n\ — 
is positive (negative) for positive (negative) total density 
n. Recalling that layer 1 is closest to the back gate, this 
shows that t he bi layer is polarised along the electric field, 
as expected [l45j . 

A single back gate in the absence of additional charged 
dopants may be described by Vt = Ubo = n to = 0, re- 
sulting in simplified expressions n = eoEbVb / (eLj,) (as in 




FIG. 7. Density-dependence of the band gap U g (n) in bi- 
layer graphene as the back-gate voltage V& is varied for a 
fixed top-gate voltage Vt [13] ■ The effective offset-density is 
no = 2[eoEtVt / (eLt) + nto]. Plots were made for screening 
parameter A = 1, using U g = \U\ji/ yL^ + Ti an< l numerical 
solution of equation (|74| . 



) and t/cxt 
simplifies 



monolayer graphene [l| 
|LT| < 7i, equation ([73j 



At high density 
ligible and U{n) 



A7in 



1 



A 



In 



n± 



(75) 



~ nj_, the logarithmic term is neg- 
A 7 in/nj_ is approximately linear in 



density. Note that the band gap U g = \U\~i\/ \fU 2 + ^\ 
tends to saturate U g — > 71, even if \U\ >• 71. At low den- 
sity, |ra| <C n±, the logarithmic term describing screening 
dominates and U g \U\ s» 2 7 i(|n|/nj_)/ ln(nj_/|n|), in- 
dependent of the screening parameter A. 

The expressions (|74I75[) for U(n) take into account 
screening due to low-energy electrons in p z orbitals using 
a simplified Hamiltonian (|67j) while neglecting oth er or- 
bitals and the effects of disorder [5jJ Il45l Il49t - |l55l ] , crys- 
talline inhomogeneity [l8| and electron-electron exchange 
and correlation. Nevertheless, there is generally good 
qualitative agreement of the dependence of U (n) on den- 
sity n predicted by equations (1741751) with density func- 
tional theory calculations 

Hill" 

and experiments in- 
cluding pho toemissio n [16l [20j and infrared spectroscopy 
[H, EiS, EH, EH- Note that the Hartree screening 
model has been gen eralised to de scribe graphene trilay- 
ers and multilayers [l47l [l56l - fl58j ■ 



IV. TRANSPORT PROPERTIES 

A. Introduction 

Bilayer graphene exhibits peculiar transport proper- 
ties due to its unusual band structure, described in sec- 
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FIG. 9. Two-probe bilayer graphene device with armchair 
edges, width W and length L. The rectangular shaded regions 
represent the ends of semi-infinite leads. 



FIG. 8. Density-dependence of the difference in layer densities 
ni — n-2 in bilayer graphene as the back-gate voltage Vb is 
varied for a fixed top-gate voltage Vt Plots were made 

for screening parameter A = 1 and the effective offset-density 
is no = 2[eo£tVt / (eLt) + nto] = 0, corresponding to the data 
labelled no = in Figure obtained by numerical solution 
of equation (|74[) . The dashed line shows the contribution 
of the partially filled bands a s me asured with respect to the 
charge neutrality point (|71[l |l48l ]. the dotted line shows the 
contribution of the f ull va lence band (|72p [la ], the solid line 
is their sum (JTJJl [l7l.[l4i|. 



tion HH where the conduction and valence bands touch 
with quadratic dispersion. Transport characteristics and 
the nature of conductivity near th e Dirac p oint were 
probed experimental ly El j,ll, [l^ . [20lll59l - [l6^ | and inves- 
tigated theoretically [l63l - ll75j | . Neglecting trigonal warp- 
ing, t he minim al conductivity is predicted to be 8e 2 /(irh) 
Il68| |. twice the value in monolayer graphene, 
while, in t he p r esen ce of trigonal warping, it is larger, 
2Ae 2 /{irh) [l63l Il69l |. because of multiple Fermi surface 
pockets at low energy, section Hi El 

For a detailed review of the electronic transport prop- 
ertie s of graphene monolayer and bilayers, see Ref. [l76l 
177| . The characteristics of bilayer graphene in the pres- 
ence of short-ranged defects and l ong-ranged charged - 
imp urities have been calculated 0, 1 1631 . Il64l . 1 1661 . Il 7ll — 
1751 ] and it is predicted that the conductivity has an 
approximately linear d epen dence on density at typical 
experimental densities [172]. At interfaces and poten- 
tial barriers, conservation of the pseudospin deg r ee o f 
freedom may influence elec tronic tra nsmission jl78l . Il 79j | , 
as in monolayer graphene [l78l Il80j | , i ncluding transmis- 



sion at monolayer-bilayer int erfac es [18lMl84| . through 
mult i ple e lectrostatic barriers [185}, or magnetic barriers 
[l86l Il87| . Inducing interlayer asymmetry and a band 
gap using an external gate [9L 1 1 7lj , describe d in sectio n lllll 
may be used to tune transport properties [l8lL Il88j | . In- 
terlayer asymmetry may also be viewed as creating an 
out-of-plane component of pseudospin and interfaces be- 



tween regions of opposite polarity have attracted theo- 
retical attention due to the existence of one- dimensio nal 
valley-polarised modes along the interface [I89l4l92j . a 
pseudospintronic analogy of spin- va lve device s for trans- 
port perpend icula r to the interface |l93l Il94| , electronic 
confi nement |l95l |. and valley- dependent transmission 

[HH. 

In the following wc describe two different models of 
the conductivity of bilayer graphene at low energy. The 
first is for ballistic transport in a clean device of finite 
length that is connected to semi-infinite leads, described 
using wave matching to calculate the transmission prob- 
ability and, then, the conductance. The second model 
describes the conductivity of a disordered, infinite sys- 
tem using the Kubo formula and the self-consistent Born 
approximation to describe scattering from the disordered 
potential. Although the two models are quite different, 
both predict t he minimal conductivity to be 8e 2 /(irh) 
[HI [163, EH- Finally, in section IT7C21 we describe 
localisation effects. 



B. Ballistic transport in a finite system 

Ballistic transport in a finite, mesoscopic bilayer 
graphene nanostructure has been modelle d in a number 
of papers flM [167^1691 li~78l EM EMEM ]. Here, we fol- 
low a wave-matching approach of Snyman and Beenakker 
|l68j |. For bilayer, as compared to monolayer, there is a 
new length scale l\ — hv/ji characteristic of the inter- 
layer coupling. Here, v — v / 3a7o/2/i is the band velocity 
of monolayer graphene, so l\ = (y/3/2) (70/71)0 « 18 A 
is several times longer [l68| than the lattice constant 
a = 2. 46 A [53| ■ For most situations, the sample size 
I >4i and the device generally behaves as a (co upled) 
bilayer rather than two separate monolayers [l68| . 

We consider a two-probe geometry with armchair edges 
as shown in Fig. |H1 There is a central mesoscopic bilayer 
region, width W and length L, connected to a left and 
right lead. This orientation is rotated by 90° as compared 
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to that described in section |TT] so the corners of the Bril- 
louin zone are located at wavevectors = £(0, 47r/3a). 
In terms of wavevector measured from the centre of the 



valley, 



ky +£47r/3a, then the Hamiltonian (j3"0)) 



in basis Al, Bl, A2, B2 may be written as 



where tt 



-ih{k x 



( U V7rt 
V7T U 

7i 

V o o 





7i 






U VTT^ 
VTT U 



\ 



(76) 



i£,k y ), 7r^ = ih(k x — i£k y ), and U is 
the on-site energy which describes the doping of the bi- 
layer. For simplicity, we include only the main interlayer 
coupling term 71 of the orbitals on the dimer Bl and A2 
sites. It is assumed that the transverse wavevector k y 
is real and it is conserved at the interfaces between the 
bilayer and the leads. The Hamiltonian (|76|) shows there 
are two values of the longitudinal wavevector k x for given 
U, k y and energy e, which we denote as fc+ and 



hvk± 



h 2 v 2 k 2 v . 



(77) 



Left- $± and right- $^ moving wave functions may be 
written as 



$^ = N± 



$f = N± 



^ ^fihv(—k± — i^ky) \ 
T(e - U) 
(e-U) 
\ —ihv(—k± + i£k y ) J 

^ Tifrv(k± — i^ky) \ 
T(e-U) 
(e-U) 
\ —ihv(k± + i^ky) ) 

[4W(e 



-ik±x-\-ik y y 



ik±x-\-ik y y 



■ (78) 



(79) 



f/)fc±]~ 1/2 for unit 



where normalisation N± 
current. 

The aim is to describe a mesoscopic bilayer region of 
finite length L connected to macroscopic leads. In order 
to mimick macroscopic, metallic contacts, there should 
be many propagating modes in the leads that overlap 
with the modes in the central bilayer region. If this is 
the case, then the value of minimal conductance should 
not depend on the particular model used for the leads, 
e.g., square lattice or gr aphe me lattice, as demonstrated 
for monolayer graphene 1971 ] . Note, this will yield quite 
different results from a model with a bilayer lead at the 
same level of doping as the central region; then, the sys- 
tem is effectively an infinite system, not a finite, meso- 
scopic conductor. 

Snyman and Beenakker [l68j j modelled the leads as 
heavily-doped bilayer graphene, gene ralis ing an approach 
developed for monolayer graphene [198]. In this way, 
there are many conducting modes present in the leads 
and it is possible to simply use matching of the bilayer 
wave functions at the interface between the central re- 
gion and the leads. In particular, the leads are mod- 
elled as bilayer graphene with on-site energy U — —Uoo 



where U^ > and U^ >• {|e|,7i}. Then, in the leads, 
k + fa fc_ « Uoo/(tiv) and wave functions ip\ e h,±, bright, ±> 
in the left (x < 0) and right (a: > L) leads may be written 
as 
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e ik y y , (80) 



a iUoo (x — L) / hv-\-ik y y 



Here, a right-moving wave with unit flux correspond- 
ing to a state with wavevector k± ps Uoa/(hv) is in- 
jected from the left lead [the first term in equation (|80|) ]. 
Subsequently, there are two left-moving waves that have 
been reflected with amplitudes r+ and , and two right- 
moving waves are transmitted to the right, with ampli- 
tudes i+ and t_. 

At the charge-neutrality point e = U = in the central 
bilayer region, the wave functions are evanescent with 
imaginary longitudinal wavevector, equation ([77| . Two 
states with k x = —i^k y have finite amplitude only on the 
Al, A2 sites and two with k x = i^k v have finite amplitude 
only on the Bl, 52 sites: 
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For each incoming mode ± from the left lead [the first 
term in equation (I80p ]. there are eight unknown ampli- 
tudes cf , cf , c 3 , cf , r + , r*, t^_, t_. Continuity of the 
wave functions at the interfaces x — and x — L between 
the central region and the leads provides eight simulta- 
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neous equations: 
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(81) 
(82) 
(83) 
(84) 

(85) 

(86) 
(87) 

(88) 



where i\ — hv/ji. Solving yields the transmission matrix 



tx t: 
t+ t: 



2/ 



2 + (L/^i) 2 + 2cosh(2fc y L) 



(L/li - 2i) cosh (k v L) - (L/£i) sinh (£,k y L) 
(L/£i) sinh (£k y L) - {L/h + 2i) cosh (k y L) 



The transmission probability [168| is determined by the 
eigenvalues T± of the product tt^ : 



T+ = 



i 



cosh 2 {£k y L =p fe c i) ' 
In 




(89) 
(90) 



The transmission coefficients T± have the same form 
as the transmission in monolayer graphene T = 
1/ cosh 2 (k y L) [l98] but shifted by the parameter k c , as 
shown in figure 1101 

The conductance G ma y be det ermined using the 
Landauer-Buttiker formula 199, 2001 



G = 



9vg s e 



Tr(ttt), 



(91) 



where the factor of g v 9s accounts for valley and spin de- 
generacy. For a short, wide sample whose width W ex- 
ceeds its length L, W 3> L, the transverse wavevector 
may be assumed to be continuous and 



G = 



g v g s e 2 W 
h 2?r 



(T+ + T-)dk y 



2g v g s e 2 W 
h ttL 



(92) 



Thus, the minimal conductivity a = GL/W = 8e 2 /(irh) 
is twice as large as in the monolayer. In a similar way, it is 
possible to determine the Fano factor of shot noise which 
take s the same value 1/3 |l68j | as in monolayer graphene 
Transmission via evanescent modes in graphene 
has been described as pseudodiffusive because the Fan o 
factor takes the same value as in a diffusive metal Il98ll. 




FIG. 10. The transmission coefficients T± of bilayer 
graphene (HHJ for L = 50£i [Hi] (solid li nes) . The trans- 
mission coefficient of monolayer graphene [l98t ] is shown in 
the centre (dashed line). 



C. Transport in disordered bilayer graphene 

1. Conductivity 

When the Fermi energy ef is much larger than the level 
broadening caused by the disorder potential, the system 
is not largely different from a conventional metal, and the 
conductivity is well described by Boltzmann transport 
theory. However, this approximation inevitably breaks 
down at the Dirac point, where even the issue of whether 
the system is metallic or insulating is nontrivial. To 
model electronic transport at the charge neutrality point, 
we need a refined approximation that properly includes 
the finite level broadening. Here, we present a conduc- 
tivity calculatio n usin g the self-consistent Born approxi- 
mation (SCBA) [l63j |. We define the Green's function as 
G(e) = (e — Ti)^ 1 . The Green's function averaged over 
the impurity configurations satisfies the Dyson's equation 

(G a , a ,(e)) =^G(°)(£) 

+ Gl°'(^E a , (11 (e)(G ai ,.( E )), (93) 
«i 

where ( ) is an average over configurations of the disorder 
potential, a is an eigenstate of the ideal Hamiltonian Hq, 
and Gq°' = (e — e Q ) _1 , with e a being the eigenenergy of 
the s tate a in Ho- In SCBA, the self-energy is given by 



201] 



Ea.a'fc) = X! ( u °,c. l U a i,a>){G ai , ai (e)). (94) 



a 1 ,a 1 



The equations (|9"3"|) and ([M]) need to be solved self- 
consistently. The conductivity is calculated using the 
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Kubo formula, 



Jv9s 



he 2 
2tt^ 



RcTr [v x (G 



R )v RA 



(G J 



v x {G R )v RR {G R )] 
(95) 



where G R = G(e + iO) and G A = G(e - iO) are retarded 
and advanced Green's functions, v x = dHo/dp x is the 
velocity operator, and g v 9s accounts for summation over 
valleys and spins. v R and v RR the velocity operators 
containing the vertex correction, defined by v RA = v x (e+ 
iO, e - iO) and v RR =v x (e + iO, e + iO) with 



v x (e,e')=v x + {UG(e)v x G(£')U). 



(96) 



In SCBA, v x should be calculated in the ladder approxi- 
mation. 

For the disorder potential, we assume a short-ranged 
potential within each valley, 



/iooo\ 

10 
10 

V o o o i / 



(97) 



and neglect intervalley scattering between K + and K_. 
This situation is realized when the length scale of the 
scattering potential is larger than atomic scale but much 
shorter than the Fermi wave length. We assume an equal 
amount of positive and negative scatterers Ui = ±u and 
a total density of scatterers per unit area ^imp • 

The SCBA formulation is applied to the low-energy 
Hamiltonian equation (|38l) . with the trigonal warping ef- 
fect due to 73 included 163]. Energy broadening due to 
the disorder potential is characterised by 



2 lmp 2-kK 

When T>sl, i-e., the disorder is strong enough to smear 
out the fine band structure of the trigonal warping, we 
can approximately solve the SCBA equation in an ana- 
lytic form. The self-energy, equation (|94p , becomes diag- 
onal with respect to index a and it is a constant, 



S(e + i0) re -iT. 



(99) 



Then, the conductivity at the Fermi energy e is written 
as 









1 + 


(£- 









l£l 

r 



4ire L 



7T z n. 2 

(100) 

The third term in the square bracket arises from the ver- 
tex correction due to the trigonal warping effect, and £l 
is given by equation (|42)l . For high energies |e| » T, a 
approximates as 



a(e) w g v c 



e 7T |e| 



(101) 




FIG. 11. (a) Calculated density of states and (b) conduc- 
tiv ity a s a function of energy for different disorder strength 
F [2021 ] . In (b), broken curves show the results of the low- 
energy two-band model, and the horizontal line indicates the 
universal conductivity g v g s e 2 /(ir 2 h). 



which increases linearly with energy. The value at zero 
energy becomes 



o-(0) = 9v9s 



n 2 h 



r J 



(102) 



In the strong disorder regime T 3> 2tt£l, the correction 
arising from trigonal warping vanishes and the c onductiv- 
ity approaches the universal value g v gse 2 / (^h) (l63lll67 1 
which is twice as large as that in monolayer graphene in 
the same approximation. In transport measurements of 
suspended bilayer graphene jl6l| . the minimum conduc- 
tivity was estimated to be about 10 _4 S, which is close to 
g v g s e 2 /{Tt 2 h). 

The 2x2 (two-band) model works well at low energy, 
but it is not expected to be valid in the strong disorder 
regime when mixing to higher energy bands is consid- 
erable. To see this, we numerically solved SCBA equa- 
tion for the original 4x4 (four-band) Hamiltonian. Fig- 
ure [TlTa) and (b) show the density of states (DOS) and 
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cond uctivity, respectively, for several disorder strengths 
[202j |. In (b) the results for the 2x2 model in equa- 
tion (jlOOp are expressed as broken curves. In (a), we 
observe that the DOS at zero energy is significantly en- 
hanced because states at high energies are shifted toward 
the Dirac point by the disorder potential. However, the 
zero-energy conductivity barely shifts from that of the 
2x2 model [equation (|102j) ]. even for strong disorder such 
as r/71 = 0.15, where the DOS at zero energy becomes 
nearly twice as large as in the 2x2 model. For higher en- 
ergy |e| > T, the conductivity increases linearly with |e| 
in qualitative agreement with equation (|100j) of the 2x2 
model, while the gradient is generally steeper. The con- 
ductivity has a small dip at the higher band edge around 
|e| ~ 71, because the frequency of electron scattering is 
strongly enhanced by the higher band states. 

The SCBA calculation was recently extended for long- 
range scatterers [203| ] . It was shown that the conductivity 
at zero energy is not universal but depends on the degree 
of disorder for scatterers wi th lo ng-range potential, sim- 
ilar to monolayer graphene 204j . 



2. Localisation effects 

The SCBA does not take account of some quantum 
corrections, such as those included explicitly in weak lo- 
calisation. In graphene, the presence of spin-like degrees 
of freedom related to sublattices and valleys, as well as 
real electronic spin itself, creates the possibility of a rich 
variety of quantu m interference behaviour. Weak local- 
isation [205l |206| | is a particularly useful probe because 
it is sensitive to elastic scattering that causes relaxation 
of the sublattice pseudospin and valley 'spin'. In the 
absence of symmetry-breaking scattering processes, con- 
servation of pseudospin in mono layer graphene tends to 
suppress backscattering [H, l207j | , and the interference of 
chiral el ectro ns would be expected to result in antilocal- 
ization [2081 ]. However, intravalley symmetry- breaking 
relaxes the pseudospin and suppresses anti-localisation, 
while intervalley s cattering tends to restore conventional 
weak localisat ion [209|-|213|, as observed experimentally 
[209l [2li]219t . Nevertheless, a nti-l o calis ation has been 
observed at high temperature [21 8l . [220] when the rel- 
ative influence of symmetry-breaking disorder is dimin- 
ished, and its pr esence ha s also been predicted at very 
low temperature [22lM223 | when spin-orbit coupling may 
influence the spin of the interfering electrons. 

In bilayer graphene, the pseudospin turns twice as 
quickly in the graphene plane as in a monolayer, no 
suppression of backscattering is expected and the quan- 
tum correction should be conventional weak localisation 
[2121 12131 T224], However, the relatively strong trigonal 
warping of the Fermi line around each valley (described 
in section III El) can suppress localis ation unless interval- 
ley scattering is sufficiently strong [224]. Experimental 
observations confirmed this picture [1591 ] , and it was pos- 
sible to determine the temperature and density depen- 



dence of relevant relaxation lengths by comparing to the 
predicted magnetoresistance formula [224(. 

Localisation has also been studied for gapped bilayer 
graphene in the presence of interlayer potential asymme- 
try U [225t | . It was shown that, as long as the disorder 
potential is long range and does not mix K± valleys, gap 
opening inevitably causes electron delocalisation some- 
where between U — and U = 00, in accordance with 
the transition of quantum valley Hall conductivity, i.e., 
the opposite Hall conductivities associated with two val- 
leys. This is an analog of quantum Hall physics but can 
be controlled purely by an external electric field without 
any use of magnetic fields. 



V. OPTICAL PROPERTIES 

The electronic structure of bilayer graphene was 
probed by spectrosc opic measure ments in zero magnetic 
" |55|,j56|, [nl M, [TH, [Hj, and also in high mag- 
The optical absorption for perpendic- 



field [16 
netic fields [7S 
ularly incident light is described by the dynamical con- 
ductivity in a electr ic field paralle l to the layers, in both 



symmetric bilayers [HI [226H228I] and in the presence of 
an interlayer-asymmetry gap |227j For sy mmetric bilay er 
graphene, this is explicitly estimated as [l64l l226M22i 



Re a xx (uj) 



9v9s e 



2 r Huj + 271 



16 h \ hw + 71 

2 



6(huj - 2|e F |) 



+ - 71) + - 71 - 2M 



hio — 271 
hu — 71 



9{huj - 2 7 i) 



-71 log 



2\e F \ +71 



7i 



S(hu> — 71) 



(103) 



where Ef is the Fermi energy and we assumed \ef\ < 71 ■ 
We label the four bands in order of descending energy 
as 1,2,3,4. The first term in equation (|103|) represents 
absorption from band 2 to 3, the second from 2 to 4 or 
from 1 to 3, the third from 1 to 4, and the fourth from 3 
to 4 or from 1 to 2. Figure[T2"lfa) shows some examples of 
calculated dynamical conductivity Re a xx (uj) with several 
values of the Fermi energy [229]. The curve for £p = 
has essentially no prominent structure except for a step- 
like increase corresponding to transitions from 2 to 4. 
With an increase in ef, a delta- function peak appears at 
huj = 71 , corresponding to allowed transitions 3 to 4. 

In a magnetic field, an optical excitation by perpen- 
dicular incident light is only allowed between the Lan- 
dau levels with n and n ± 1 for arbitrary combinations 
of /i = H, L and s = ±1, since the matrix element of 
the velocity operator v x vanishes otherwise. Figure [T^] 
(b) shows some plots of Re a xx (uj) in magnetic fields at 
Ef = and zero temperature [230| . Dotted lines pene- 
trating panels represent the transition energies between 
several specific Landau levels as a continuous function 
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FIG. 12. Interband part of the dynamical conductivity of 
bilayer graphene plotted against t he fr equency cj, in (a) zero 
magnetic field with diffe rent £f's [229] (b) several magnetic 
fields with ef = [2301 ]. Dashed curves in (b) indicate the 
transition energies between several Landau levels. 



of TujJb- Every peak position behaves as a linear func- 
tion of B (x hujg in weak field but it switches over to 
V^B-dcpendence as the corresponding energy moves out 
of the parabolic band region. In small magnetic fields, 
the peak structure is smeared out into the zero-field curve 
more easily in the bilayer than in the monolayer, because 
the Landau level spacing is narrower in bilayer due to the 
finite band mass. 



VI. ORBITAL MAGNETISM 

Graphite and related materials exhibit a strong or- 
bital diamagnetism which overcomes the Pauli spin para- 
magnetism. Theoretically the diamagnetic s usceptibil- 
ity was calculated for gra phite Il2ll l23ll |232| , graphite 
intercalation compounds as well as few-layer 



graphenes [237, |238|. In particular, monolayer graphene 
has a strong diamagnetic singularity at Dirac point, 
whi ch is expresse d as a Delta function in Fermi energy 
£f [l2l[ |239| - |242| ] . In the bilayer, the singularity is re- 
laxed by the modification of the band structure caused 
by the in terlayer coupling as we will see in the following 
f234Ll237t 

When the spectrum is composed of discrete Landau 
levels £„, the thermodynamical potential is generally 
written as 



1 



9v9s 



/3 2tt1 2 



E 



(104) 



where /3 = l/k B T, ip(e) = log [l + e- 13 ^"^] with ( being 
the chemical potential. In weak magnetic field, using the 
Euler-Maclaurin formula, the summation in n in equa- 
tion (|104[) can be written as an integral in a continuous 
variable x with a residual term proportional to B 2 . The 
magnetization M and the magnetic susceptibility x can 
be calculated by 



M 



sdn\ 

AdBJc' 



X 



dM 
~dB 



B=0 



/d 2 Q\ 



B=0 



(105) 



For monolayer graphene, the susceptibility is [121LI23 

dm 



X 



-9v9s 



2 2 
e v 

6nc 2 



de 



de. 



(106) 



At zero temperature, it becomes a delta function in Fermi 
energy, 



x(£f) = -g v gs 



6ttc 2 



S(e F ). 



(107) 



The delta-function susceptibility of monolayer graphene 
is strongly distorted by the interlayer coupling 71. For 
the Hamiltonian of the symmetric bi layer grap hene, the 
orbital susceptibility is calculated as [234 l237j ] 



X(e) = 9v9s 



2 2 
e v 

47TC 2 7] 



■0(7i 



N)(iogM + I 



(108) 



The susceptibility diverges logarithmically at ep — 0, be- 
comes slightly paramagnetic near \ef\ = 7i, and vanishes 
for \ef\ > 71 where the higher subband enters. The in- 
tegration of x i n equation (|108j) over the Fermi energy 
becomes — g v g s e 2 v 2 / (3irc 2 ) independent of 71, which is 
exactly twice as large as that of the monolayer graphene, 
equation (|107[) . 

The susceptibility wa s als o calculated in the presence 
of interlayer asymmetry [2431 ] . Figure[13](a) and (b) show 
the density of states and the susceptibility, respectively, 
for bilayer graphene with {7 = 0, 0.2, and 0.5. The 
susceptibility diverges in the paramagnetic direction at 
the band edges where the density of states also diverges. 
This huge paramagnetism can be interpreted as the Pauli 
paramagnetism induced by the valley pse udo-spin split- 
ting and diverging density of states [244|. The suscep- 
tibility vanishes in the energy region where the higher 
subband enters. 
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FIG. 13. (a) Density of states, and (b) susceptibility of bilayer 
grap henes with the asymmetry gap U/ji = 0, 0.2, and 0.5 
[243l ]. In (b), the upward direction represents negative (i.e., 
diamagnetic) susceptibility. 



VII. PHONONS AND STRAIN 

A. The influence of strain on electrons in bilayer 
graphene 

Deformation of a graphene sheet couples to the elec- 
tronic system and modifies the low-energy Hamiltonian. 
In monolayer graphene, static changes in distance and 
angles of the atomic bonds can be described as effec- 
ti ve sca lar or vector potentials in the Dirac Hamiltonian 
@; l245j . Bilayer graphene has extra degrees of freedom in 
deformation associated with the presence of two layers. 
It was shown that a band gap can be opened by giving 
different distortions to the two layers or by p ulling th e 
two layers apart in the perpendicular direction [246l - l248j . 

Rather than produce a band gap, it has been predicted 
that homogeneous lateral strain in bilayer graphene can 
produc e a change in topology of the band structure at low 
249M251) . formation causes tight-binding 



energy 

parameters 70 and 73 to become dependent on the hop- 
ping direction and produces an additional term in the 
low-energy two-band Hamiltonian (|38|) 



(109) 



where parameter w = \w\e % ^ 6 depends on the microscopic 




details of the deformation; it is non-zero only when the 
role of skew interlayer coupling 73 is taken into account 
[249j . It is estimated that 1% strain would give \w\ ~ 
6meV [249( 1 - Taken with the quadratic term ho m the 
Hamiltonian (|38p , the low-energy bands with energy E — 
±£i are given by 



\w\p- 



cos (2tp + 9) 



h) ■ (110 » 



This describes a Lifshitz transition at energy sl = \w\, 
below which there are two Dirac point s in the vicinity 

1914251 



24C 



centred at mo- 
\/\w\ 71/w and angles ip = —9/2 



of each Brillouin zone corner 
mentum p sa \j2m \w\ 
and ip = tt — 9/2. In general, there should be an interplay 
between terms h s , equation (|109j) . and h w , equation (|38|) . 
leading to the possibility of employing strain to annihi- 
late two Dirac points and, thus, change the low e nerg 
topology of the bands from four to two Dirac points (24S 

The presence of two Dirac points would cause zero- 
energy Landau levels to be eightfold degenerate; an ex- 
perimental signature of this state is predicted to be the 
persistence of fill ing f actor v = ±4 in the low-field quan- 
tum Hall effect [2491 ] . This contrasts with the Lifshitz 
transition that would occur in the presence of parameter 
73 without strain when there are four Dirac points, sec- 
tion III El giving a degeneracy of sixteen and v = ±8 at 
low fields. In both cases, Berry phase 27r is co nserved: 
two Dirac cones with Berry phase 7r eac h 12491 or four 
Dirac cones with three of 7r and one of — 7r [12|,|94|]. ^ has 
also been predicted that the presence of the Lifshitz tran- 
sition will be notice able in the low-energy conductivity at 
zero magnetic field [l63l Il69| | , and the particular case of 
two Dirac points in the presence of strain has recently 
been analysed, too [252|. Note that the effect of lateral 
strain on the low-energy topology of the band structure 
is qualitatively similar to that of a gapless nematic phase 
which possibly arises as the re sult of el ectron-electron 
interactions in bilayer graphene 253j-|256|. 



B. Phonons in bilayer graphene 

Raman spectroscopy has been a valua ble t ool in prob- 
ing the behaviour of phonons in graphite [2571 ] and it may 
be used t o det ermine the number of layers in multilayer 
graphene 258], differentiating between monolayer and bi- 
layer. For an in-depth review of Raman spectroscopy 
of graphe ne includ ing bilayer graphene see, for exam- 
ple, Refs. 259, 260]. The phonon spectrum of monolayer 



graphene has been calculated using a tight-binding force- 
constant model with parameters fit to Raman data [26 lj ] , 
and with density functional theory [262h265| . There are 
three acoustic (A) and three optical (O) branches consist- 
ing of longitudinal (L) and transverse (T) in-plane modes 
as well as out-of-plane (Z) modes. At the zone centre (the 
r point), the TA and LA modes display linear dispersion 
uj <~*j q but the ZA mode is quadratic u> ~ q 2 . The ZO 
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mode is at ~ 890 cm [26 



34], and the LO and TO 
modes are degenerate (at ~ 1600 cm -1 ). At the K point, 
the ZA and ZO modes (~ 540 cm -1 ) and the LA and LO 
modes (~ 1240 cm -1 ) are degenerate, with TA modes 
at ~ 1000 cm -1 and TO at - 1300 cm -1 . For undoped 
graphene, owing to strong electron-phonon coupling, the 
highest optical modes at the T and K point (i.e. the LO 
mode at the F point and the TO mode at the K point) 
display Kohn anomalies [26614268} whereby the phonon 
dispersion to(q) has an almost linear slope as observed, 
for examp le, in inelastic x-ray measurements of graphite 
[263l l268j |. As graphene is a unique system in which the 



electron or hole concentration can be tuned by an ex- 
ternal gate voltage, it was realised [269L l27Cl| that the 
change in electron density would also influence the be- 
haviour of the optical phonons through electron-phonon 
coupling and, in particular, a logarithmic singularity 
in their dispersion was predicted |269j when the Fermi 
energy sp is half of the energy of the optical phonon 
|ejr| = huj/2. Subsequently, such tuning of phonon fre- 
quency and bandwidth by adjusting the electronic den- 
sity was obser ved i n mo nolayer graphene through Raman 
spectroscopy [27l[ |272| . 

The behaviour of phonons in bilayer graphene has been 
observe d exp erimentally through Raman spectro scopy 
[2i [z3, f273M279| j and infrared spectroscopy [US [28lf . 
with particular focus on optical phonon anomalies and 
the influence of gating. Generally, the phonon spectrum 
of bilayer graphene, w hich has b een calculated using den- 
sity functional theory [2821 1283| and force-constant mod- 
els [284 |285| . is similar to that of monolayer. Near the 
r point there are additional low-frequency modes. There 
is a doub ly-degenerate rigid shear mode at ~ 30 c m -1 
[284 |285| observed through Raman spectroscopy [278j | 
and an optical mode at ~ 90 cm -1 which arises from 
relative motion of the layers in the vertical direction 
(perpendicular to th e layer plane), known as a layer- 
breathin g^ mo de [286| and observed through Raman spec- 
troscopy [279j . At the T point, inter layer coupling causes 
the LO/TO modes to split into two doubly-degenerate 
branches where the higher (lower) frequency branch cor- 
responds to symmetric 'in-phase' (antisymmetric 'out-of- 
phase') relative motion of atoms on the two layers (in the 
in-plane direction). Analogously to the monolayer, it was 
predicted that these optical phonons would be affected 
by electron-phonon coupling, with a logarithmic singu- 
larity in the dispersion of the symmetric modes when the 
Fermi ene rgy g p is equal to half of the optical phonon 
frequency 287], and hybridisation of the symmetric and 
antisymmetric mode s in t he presence of interlayer poten- 
tial asymmetry [288l l289j | . Experimentally, this anoma- 
lous phonon dispersio n has bee n observed through Ra- 
man spectroscopy [77], I273I4276I ] including the evolution 
of two distinct components i n the Ram an G band for non- 
zero interlayer asymmetry [274 - 1276 1. The Raman spec- 
trum has also been studied for bilayer grap hene in th e 
presence of Landau levels in a magnetic field [287, 290]. 



C. Optical phonon anomaly 

In the following, we describe the anomalous optical 
phonon spectrum in bilay er graphe ne taking into account 
electron-phonon coupling [287l . l288l |. Theoretically, it was 
shown that a continuum model works wel l in describ- 
ing long- wave length acoustic phonons [245| and optical 
phonons [29 ll ] in gr aphene, an d this theory was extended 
to bilayer graphene [287l l288j | . An optical phonon on one 
graphene layer is represented by the relative displacement 
of two sub-lattice atoms A and B as 



u(r) = 



ANMuq 



(& q ,M + &L,Je M (q)e^' r , (111) 



where N is the number of unit cells, M is the mass of a 
carbon atom, u>o is the phonon frequency at the V point, 
Q = ilx^Qy) is the wave vector, and 6 qjA1 and &], are 
the creation and destruction operators, respectively. The 
index \x represents the modes (t for transverse and I for 
longitudinal) , and corresponding unit vectors are defined 
by ei = iq/|q| and e t = iz x q/|q|. 

The Hamiltonian of optical phonons is written as 



P h 



(112) 



and the interaction with an electron at the K + point is 



H 



K 
i nt 



(3hv 

[a x Uy(v) - a y u x (r)] 

a cc 



(113) 



where the Pauli matrix cr, works on the space of 
(4>A\, 4>bi) for the phonon on layer 1, and (<pA2, 4>B2) for 
layer 2. The dimensionless parameter /3 is related to the 
dependence of the hopping integral on the interatomic 
distance, and is defined by (3 = — dlogj /dlogacc- We 
usually expect (3 ~ 2. The strength of the electron- 
phonon interaction is characterized by a dimensionless 
parameter 



36^3 h 1 f/3 



2Mb 2 hion 



(114) 



For M = 1.993 x 10 23 g and haj = 0.196eV (correspond- 
ing to 1583cm -1 ), we have A 3 x 10 -3 (^/2) 2 . For 
the K- point, the interaction Hamiltonian is obtained 
by replacing with — a\. 

The Green's function of an optical phonon is given by 
a 2 x 2 matrix associated with phonons on layers 1 and 
2. This is written as 



D(q,w) 



(Huj) 2 - (huj ) 2 - 2fiw n(q,a;) 
1 

hui — huo — n(q, luq) 



(115) 



where il(q, uS) is the phonon's self-energy, and the near- 
equality in the second line stands because the self-energy 
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FIG. 14. (a) Frequency shift and (b) broa denin g of the T- 
point optical phonon in a bilayer graphene 288]. Solid and 
dashed lines denote the high- and low- frequency modes, re- 
spectively, and the thin dotted lines in the top panel show 
the frequencies for symmetric and antisymmetric modes cal- 
culated without inclusion of their mixing. 



is much smaller than Hloq. Then, the eigenmodes are 
given by eigenvectors \u) of ReII(q, cJrj), and the fre- 
quency shift and broadening are obtained as the real and 
imaginary part of (it|II(q, uq)\u). 

In symmetric bilayer graphene (i.e., no inter layer po- 
tential difference), eigenmodes are always classified into 
symmetric and antisymmetric modes in which the dis- 
placement of the top and bottom layers is given by (u, u) 
for the former, and (u, — u) for the latter. The symmet- 
ric mode causes interband transitions between the con- 
duction band and valence band, while the antisymmetric 
mode contributes to the transitions between two conduc- 
tion bands (or between two valence bands). The potential 
difference between two layers gives rise to hy bridisation 
of symmetric and antisymmetric modes 287]. 

Figure [14] (a) and (b) show the calculated frequency 
shift and broad ening , respectively, as a function of elec- 
tron density n s [2881 ]. Here we take hajo = 71 /2, and in- 
troduce a phcnomenological broadening factor of O.lfiwrj. 
We assume that n s is changed by the gate voltage, and 
appropriately take account of the band deformation due 
to the interlayer potential difference when n s 7^ 0, with 
the self-consistent screening effect included. The thin 
dotted lines in the top panel in Figure [14] (a) indicate 
the frequencies for symmetric and antisymmetric modes 
calculated without inclusion of their mixing. The lower 
and higher branches exactly coincide with symmetric and 
asymmetric modes at n s = where the mixing is absent. 
The dip at the symmetric mode occurs when huo = 2e^, 



i.e., the interband transition excites a valence electron 
exactly to the Fermi surface. The coupling between sym- 
metric and antisymmetric modes arises when n s ^ 0, and 
makes an anti-crossing at the intersection. 



VIII. ELECTRONIC INTERACTIONS 

Generally speaking, the low-energy behaviour of elec- 
trons in bilayer graphene is well described by the tight- 
binding model without the need to explicitly incorporate 
electron-electron interaction effects. Coulomb screening 
and collective excitations have been des cribed in a num- 
ber of theoretical papers [292ll300t and the impor- 
tance of interaction effects in a bilayer under external 
gating [lil [30lM307t has been stressed. Interaction ef- 
fects should also be important in the presence of a mag- 
netic field or at very low carrier density, particularly in 
clean samples. 

Bilayer graphene has quadratic bands which touch at 
low energy resulting in a non-vanishing density of state 
and it has been predicted to be unstable to electron- 
electron interactions at half filling. Trigonal warping 
tends to cut off infrared singularities and, thus, finite 
coupling strength is generally required to realise corre- 
lated ground states; if trigonal warping is neglected, then 
arbitrarily weak interactions are sufficient. Since bilayer 
graphene possesses pseudospin (i.e., which layer) and val- 
ley degrees of freedom, in addition to real electron spin, it 
is possible to imagine a number of different broken sym- 
metry states that could prevail depending on model de- 
tails and para meter values. Suggestions inclu de a ferro - 
magnetic I308J . la yer an tiferromagnetic [ill 1^ . 1309143 12 1. 
ferro electric |31ll . |313| or a charge density wave state 
[314j ] ; topologically non-trivial phases with bulk gaps and 
gapless ed ge excitations suc h as an anomalous quantum 
Hal l stat e [87ll31ll . [315ll316l ] or a quan tum spin Hall state 
[87t 1 3 (also called a spin fl ux state |256| ); or a gapless 
nematic phase [253142561 l317j | . 

Insulating states contribute a term proportional to a z 
in the two-band Hamiltonian (|38|) with its sign corre- 
sponding to the distribution of layer, valley and spin 
degrees of freedom, indicated in figure [T5] as manifest 
in their spin - and valley-dependent Hall conductivities 
[13, H3, l31l| . Note that the quantum spin Hall state 
produces a term equivalent to that of intrinsic spin-orbit 
coupling, equatio n (HoT). w hich describes a topological in- 
sulator state (99l . lllOl . il 111 ]. By way of contrast, the gap- 
less nematic phase has a qualitatively similar effe ct on 
the spectrum of bilayer graphene as lateral strain [249( , 
described in section |VI"TI producing an additional term in 
the two-band Hamiltonian of the form of equation (|109[) . 
with parameter w taking the role of an order parameter. 

Experiments on suspended bilayer graphene devices 
have found evidence for correla ted state s at very low 
density and zero magnetic field 31814324 1 . Conductiv- 
ity measurements of double-gate devices 111 observed a 
non-monotonic dependence of the resistance on electric 
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FIG. 15. Schematic of the distribution of spin, valley and 
layer degrees of freedom in candid ates for spontaneous gapped 
states in bilayer graphene [TJ, |87| . l31ll | . They each contribute 
a term proportional to cr z in the two-band Hamiltonian (|38|l , 
with its sign depending on the valley £ = ±1 and spin s = 
±1 orientation. Arrows indicate spin orientation, located at 
valley K+ (solid) or K- (dashed) and on different layers (top 
or bottom) of a state at the top of the valence band. 



field cumulating in a non-divergent resistance at zero field 
while compressibility measurements of single-gate devices 
|319j | found an incompressible region near the charge neu- 
trality point. These measurements were interpreted as 
being consistent with either the anomalous quantum Hall 
state, in which the separate layers of the device are valley 
polarised, or the gapless nematic phase. 

Subsequen tly, c onductivity measurements on single- 
gate devices [320l | observed a weak temperature depen- 
dence of the width of the conductivity minimum near 
zero carrier density, suggesting a suppressed density of 
states as compared to that expected for a parabolic dis- 
persion, as well as particularly robust cyclotron gaps at 
filling factor v = ±4, observations attributed to the pres- 
ence of the nematic phase. 

However, other experiments 32114324 observed insu- 
lating states indicating the formation of ordered phases 
with energy gaps of about 2 meV. There was evidence 
for the exis tence of edge states in one of th ese exper- 
iment [32l| . but not in the others and the 
latter observations i ncluding their response to perpen- 
dicular electric field [322L |323| and tilted magnetic field 
|324j j were seen as being consistent with a layer antiferro- 
magnetic state, in which the separate layers of the device 
are spin polarised, figure IT57 d). At present, there appear 
to be contradicting experimental and theoretical results, 
but it should be n oted that renormalisation group stud- 
ies |255l 12561 |325| have highlighted the sensitivity of the 



phase diagram of bilayer graphene to microscopic details. 

In the absence of interactions, the low-energy Landau 
level spectrum of bilayer graphene [i| consists of a series 
of fourfold (spin and valley) degenerate levels with an 
eightfold-degenerate level at zero energy, as described in 



section III Jl The resulting Hall conductivity consists of 
a series of plateaus at conventional integer positions of 
4e 2 //i, but with a double-sized step of 8e 2 /h across zero 
density Q. Interaction effects are expected to lift the 
level degeneracy of qu antum H all ferromagnet states at 
integer filling factors Indeed, an insulating 

state at filling factor v — and complete splitting of 
the eightfold-degenerate level at zero energy have been 
observed with quantum states at filling factors v = 0, 1, 2 
and 3 in high-mobility suspended bilayer gr aphene a t low 
fields (with all states resolved at B = 3 T) [l6ll l332j | and 
in samples on silicon substrates at high fields, typically 
above 20 T [3311331. 



The fractional quantum Hall effect has been observed 
rece ntly in mo nolayer graphene, both in suspe nded sam- 
ples [335], HH and graphene on boron nitride 337], and 
there is evidence for it in bilayer graphene 338], too. 
Strongly-correlated states at fractional filling factors and 
the prospect of tuning their prop erties has been the fo- 
cus of recent theoretical attention [339143431 ]. Clearly, the 
nature of the electronic properties of bilayer graphene in 
high-mobility samples is a complicated problem, and it is 
likely to be an area of further intense experimental and 
theoretical investigation in the following years. 



IX. SUMMARY 

This review focused on the single-particle theory 
of electrons in bilayer graphene, in the shape of the 
tight-binding model and the related low-energy effective 
Hamiltonian. Bilayer graphene has two unique proper- 
ties: massive chiral quasiparticles in two parabolic bands 
which touch at zero energy, and the possibility to con- 
trol an infrared gap between these low-energy bands by 
applying an external gate potential. These features have 
a dramatic impact on many physical properties of bi- 
layer graphene including some described here: optical 
and transport properties, orbital magnetism, phonons 
and strain. A number of topics were not covered here in 
great detail or at all; we refer the reader to relevant de- 
tailed reviews of graphene including electronic transport 
[1761 Il77j |. electronic and p hoton ic devices |344j . scan- 
ning tunn elling micro scopy [345j . Raman spectroscopy 
[259l l260ll . m agnetism |346j . spintronics and pseudospin- 



tronics [3471 ] . Andreev reflection at th e int erface with a 
superconductor and Klein tunnelling [348| . growth and 
applica tion s [2l| , and the properties of graphene in gen- 
eral @, l349j ] . Finally, although the central features of the 
single-particle theory are already established, the same 
can not be said of the influence of electronic interactions, 
which is likely to remain a subject of intense research, 
both theoretical and experimental, in the near future at 
least. 

The authors gratefully acknowledge colleagues for 
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T. Ando and V. I. Fal'ko, and we acknowledge funding 
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